Load data
Check samples distribution
dim(metadata)## [1] 50 30
library(ggplot2)
library(ggsci)
# plot
samples_plot <- ggplot(metadata, aes(x = ID, y = ..count.., fill = Depth_cm)) +
geom_bar(position = "dodge") +
facet_wrap(~ Site, scales = "free_x") +
labs(title = "Samples by Site and Depth", y = "Number of Samples", x = "Site") +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, hjust = 1)) + # Rotate labels
scale_fill_lancet()
#show
samples_plot## Warning: The dot-dot notation (`..count..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(count)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Create Phyloseq
library(phyloseq)
library(qiime2R)
#create phyloseq object
ps <- qza_to_phyloseq(
features = "../results/04.qiime/ASV_table_filter_freq218_emc.qza",
#tree = "../results/04.qiime/rooted-tree-iqtree.qza",
taxonomy = "../results/04.qiime/taxonomy.qza",
metadata = "../data/metadata.tsv")ps## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 7097 taxa and 50 samples ]
## sample_data() Sample Data: [ 50 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 7097 taxa by 7 taxonomic ranks ]
library(devtools)## Loading required package: usethis
library(microbiome)##
## microbiome R package (microbiome.github.com)
##
##
##
## Copyright (C) 2011-2022 Leo Lahti,
## Sudarshan Shetty et al. <microbiome.github.io>
##
## Attaching package: 'microbiome'
## The following object is masked from 'package:ggplot2':
##
## alpha
## The following object is masked from 'package:base':
##
## transform
library(microbiomeutilities)## Loading required package: dplyr
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
## Warning: replacing previous import 'ggplot2::alpha' by 'microbiome::alpha' when
## loading 'microbiomeutilities'
##
## Attaching package: 'microbiomeutilities'
## The following object is masked from 'package:microbiome':
##
## add_refseq
#check data
microbiome::summarize_phyloseq(ps)## Compositional = NO2
## 1] Min. number of reads = 664132] Max. number of reads = 6748463] Total number of reads = 90064984] Average number of reads = 180129.965] Median number of reads = 171486.57] Sparsity = 0.5668845991263916] Any OTU sum to 1 or less? NO8] Number of singletons = 09] Percent of OTUs that are singletons
## (i.e. exactly one read detected across all samples)010] Number of sample variables are: 29SampleIDDepthDepth_cmSiteSample_typeDateLocationTemp...C..intSalnidad..ppm.intSalinidad..ppt.intpH.intpHsupRedox.medido.en.campo1Redox.mV..Eh.intTemperatura...C.supSalnidad..ppm.supSalinidad..ppt.supRedox.medido.en.campoRedox.mV..Eh.supS.2..mg.l..campoml.de.muestra1Factor.de.dilucion1S.2..mg.l.ID_MSO.....mg.l..campoml.de.muestra2Factor.de.dilucion2SO.....mg.l.2
## [[1]]
## [1] "1] Min. number of reads = 66413"
##
## [[2]]
## [1] "2] Max. number of reads = 674846"
##
## [[3]]
## [1] "3] Total number of reads = 9006498"
##
## [[4]]
## [1] "4] Average number of reads = 180129.96"
##
## [[5]]
## [1] "5] Median number of reads = 171486.5"
##
## [[6]]
## [1] "7] Sparsity = 0.566884599126391"
##
## [[7]]
## [1] "6] Any OTU sum to 1 or less? NO"
##
## [[8]]
## [1] "8] Number of singletons = 0"
##
## [[9]]
## [1] "9] Percent of OTUs that are singletons \n (i.e. exactly one read detected across all samples)0"
##
## [[10]]
## [1] "10] Number of sample variables are: 29"
##
## [[11]]
## [1] "Sample" "ID" "Depth"
## [4] "Depth_cm" "Site" "Sample_type"
## [7] "Date" "Location" "Temp...C..int"
## [10] "Salnidad..ppm.int" "Salinidad..ppt.int" "pH.int"
## [13] "pHsup" "Redox.medido.en.campo1" "Redox.mV..Eh.int"
## [16] "Temperatura...C.sup" "Salnidad..ppm.sup" "Salinidad..ppt.sup"
## [19] "Redox.medido.en.campo" "Redox.mV..Eh.sup" "S.2..mg.l..campo"
## [22] "ml.de.muestra1" "Factor.de.dilucion1" "S.2..mg.l."
## [25] "ID_M" "SO.....mg.l..campo" "ml.de.muestra2"
## [28] "Factor.de.dilucion2" "SO.....mg.l."
Search Phylum with low prevalence
# Explore prevalence
## Get prevalence
prevdf = apply(X = otu_table(ps),
MARGIN = ifelse(taxa_are_rows(ps), yes = 1, no = 2),
FUN = function(x){sum(x > 0)})
## Add taxonomy
prevdf = data.frame(Prevalence = prevdf,
TotalAbundance = taxa_sums(ps),
tax_table(ps))
## Check prevalence at Phylum level
plyr::ddply(prevdf, "Phylum", function(df1){cbind(mean(df1$Prevalence),sum(df1$Prevalence))})plyr::ddply(prevdf, "Genus", function(df1){cbind(mean(df1$Prevalence),sum(df1$Prevalence))})## Get genus prevalence plot
prevalence_genus = subset(prevdf, Genus %in% get_taxa_unique(ps, "Genus"))
prev_genus <- ggplot(prevalence_genus, aes(TotalAbundance,
Prevalence /nsamples(ps),color=Genus)) +
geom_hline(yintercept = 0.05, alpha = 0.5, linetype = 2) +
geom_point(size = 1.5, alpha = 0.7) + theme_bw() +
scale_x_log10() + xlab("Total Abundance") + ylab("Prevalence") +
facet_wrap(~Phylum) + theme(legend.position="none")
prev_genusClean names
#check names
head(phyloseq::tax_table(ps))## Taxonomy Table: [6 taxa by 7 taxonomic ranks]:
## Kingdom Phylum Class
## ASV_2303 "d__Bacteria" "Proteobacteria" "Alphaproteobacteria"
## ASV_6775 "d__Bacteria" "Proteobacteria" "Gammaproteobacteria"
## ASV_2543 "d__Bacteria" "Chloroflexi" "Dehalococcoidia"
## ASV_4404 "d__Bacteria" "Chloroflexi" "Anaerolineae"
## ASV_1529 "d__Bacteria" "Latescibacterota" "Latescibacterota"
## ASV_779 "d__Bacteria" "Chloroflexi" "Anaerolineae"
## Order Family Genus
## ASV_2303 "Rhizobiales" "Xanthobacteraceae" "Pseudolabrys"
## ASV_6775 "Burkholderiales" "Sutterellaceae" "uncultured"
## ASV_2543 "MSBL5" "MSBL5" "MSBL5"
## ASV_4404 "SBR1031" "SBR1031" "SBR1031"
## ASV_1529 "Latescibacterota" "Latescibacterota" "Latescibacterota"
## ASV_779 "Anaerolineales" "Anaerolineaceae" "uncultured"
## Species
## ASV_2303 NA
## ASV_6775 "uncultured_bacterium"
## ASV_2543 "uncultured_bacterium"
## ASV_4404 "uncultured_bacterium"
## ASV_1529 NA
## ASV_779 "uncultured_Chloroflexi"
tail(phyloseq::tax_table(ps))## Taxonomy Table: [6 taxa by 7 taxonomic ranks]:
## Kingdom Phylum Class
## ASV_4505 "d__Bacteria" "MBNT15" "MBNT15"
## ASV_5666 "d__Bacteria" "Chloroflexi" "Dehalococcoidia"
## ASV_108 "d__Bacteria" "Nitrospirota" "Thermodesulfovibrionia"
## ASV_3882 "d__Bacteria" "Firmicutes" "Bacilli"
## ASV_6154 "d__Bacteria" "Methylomirabilota" "Methylomirabilia"
## ASV_1246 "d__Bacteria" "Proteobacteria" "Gammaproteobacteria"
## Order Family Genus
## ASV_4505 "MBNT15" "MBNT15" "MBNT15"
## ASV_5666 "t0.6.f" "t0.6.f" "t0.6.f"
## ASV_108 "uncultured" "uncultured" "uncultured"
## ASV_3882 "Lactobacillales" "Enterococcaceae" "Enterococcus"
## ASV_6154 "Rokubacteriales" "Rokubacteriales" "Rokubacteriales"
## ASV_1246 "Burkholderiales" "Comamonadaceae" "Delftia"
## Species
## ASV_4505 "uncultured_prokaryote"
## ASV_5666 "uncultured_Chloroflexi"
## ASV_108 "uncultured_Nitrospirae"
## ASV_3882 NA
## ASV_6154 "uncultured_bacterium"
## ASV_1246 NA
Agglomerate levels
#Aggregate at genus level.
ps_family <- phyloseq::tax_glom(ps, "Family", NArm = TRUE)
ps_genus <- phyloseq::tax_glom(ps, "Genus", NArm = TRUE)
ps_species <- phyloseq::tax_glom(ps, "Species", NArm = TRUE)ps_genus <- filter_taxa(ps_genus, function(x) sum(x) > 1, TRUE)
head(phyloseq::tax_table(ps_genus))## Taxonomy Table: [6 taxa by 7 taxonomic ranks]:
## Kingdom Phylum Class
## ASV_3181 "d__Bacteria" "Myxococcota" "Polyangia"
## ASV_4916 "d__Archaea" "Nanoarchaeota" "Nanoarchaeia"
## ASV_1765 "d__Bacteria" "Planctomycetota" "Phycisphaerae"
## ASV_1959 "d__Bacteria" "Proteobacteria" "Gammaproteobacteria"
## ASV_7171 "d__Bacteria" "Armatimonadota" "Fimbriimonadia"
## ASV_2739 "d__Bacteria" "Verrucomicrobiota" "Verrucomicrobiae"
## Order Family
## ASV_3181 "VHS-B3-70" "VHS-B3-70"
## ASV_4916 "Woesearchaeales" "Woesearchaeales"
## ASV_1765 "Phycisphaerales" "Phycisphaeraceae"
## ASV_1959 "Gammaproteobacteria_Incertae_Sedis" "Unknown_Family"
## ASV_7171 "Fimbriimonadales" "Fimbriimonadaceae"
## ASV_2739 "Opitutales" "Opitutaceae"
## Genus Species
## ASV_3181 "VHS-B3-70" NA
## ASV_4916 "Woesearchaeales" NA
## ASV_1765 "AKYG587" NA
## ASV_1959 "uncultured" NA
## ASV_7171 "Fimbriimonadaceae" NA
## ASV_2739 "Opitutus" NA
Change id by tax name
#Substitute these IDs with names of genus
ps_genus <- format_to_besthit(ps_genus) #no me gusta el formato, pero sirvehead(phyloseq::tax_table(ps_genus))## Taxonomy Table: [6 taxa by 8 taxonomic ranks]:
## Domain Phylum
## ASV_3181:g__VHS-B3-70 "d__Bacteria" "Myxococcota"
## ASV_4916:g__Woesearchaeales "d__Archaea" "Nanoarchaeota"
## ASV_1765:g__AKYG587 "d__Bacteria" "Planctomycetota"
## ASV_1959:g__uncultured "d__Bacteria" "Proteobacteria"
## ASV_7171:g__Fimbriimonadaceae "d__Bacteria" "Armatimonadota"
## ASV_2739:g__Opitutus "d__Bacteria" "Verrucomicrobiota"
## Class
## ASV_3181:g__VHS-B3-70 "Polyangia"
## ASV_4916:g__Woesearchaeales "Nanoarchaeia"
## ASV_1765:g__AKYG587 "Phycisphaerae"
## ASV_1959:g__uncultured "Gammaproteobacteria"
## ASV_7171:g__Fimbriimonadaceae "Fimbriimonadia"
## ASV_2739:g__Opitutus "Verrucomicrobiae"
## Order
## ASV_3181:g__VHS-B3-70 "VHS-B3-70"
## ASV_4916:g__Woesearchaeales "Woesearchaeales"
## ASV_1765:g__AKYG587 "Phycisphaerales"
## ASV_1959:g__uncultured "Gammaproteobacteria_Incertae_Sedis"
## ASV_7171:g__Fimbriimonadaceae "Fimbriimonadales"
## ASV_2739:g__Opitutus "Opitutales"
## Family Genus
## ASV_3181:g__VHS-B3-70 "VHS-B3-70" "VHS-B3-70"
## ASV_4916:g__Woesearchaeales "Woesearchaeales" "Woesearchaeales"
## ASV_1765:g__AKYG587 "Phycisphaeraceae" "AKYG587"
## ASV_1959:g__uncultured "Unknown_Family" "uncultured"
## ASV_7171:g__Fimbriimonadaceae "Fimbriimonadaceae" "Fimbriimonadaceae"
## ASV_2739:g__Opitutus "Opitutaceae" "Opitutus"
## Species
## ASV_3181:g__VHS-B3-70 "g__VHS-B3-70"
## ASV_4916:g__Woesearchaeales "g__Woesearchaeales"
## ASV_1765:g__AKYG587 "g__AKYG587"
## ASV_1959:g__uncultured "g__uncultured"
## ASV_7171:g__Fimbriimonadaceae "g__Fimbriimonadaceae"
## ASV_2739:g__Opitutus "g__Opitutus"
## best_hit
## ASV_3181:g__VHS-B3-70 "ASV_3181:g__VHS-B3-70"
## ASV_4916:g__Woesearchaeales "ASV_4916:g__Woesearchaeales"
## ASV_1765:g__AKYG587 "ASV_1765:g__AKYG587"
## ASV_1959:g__uncultured "ASV_1959:g__uncultured"
## ASV_7171:g__Fimbriimonadaceae "ASV_7171:g__Fimbriimonadaceae"
## ASV_2739:g__Opitutus "ASV_2739:g__Opitutus"
Clean rare characters
# Remove characters
taxa_names(ps_genus) <- gsub(" ", ".", taxa_names(ps_genus))
taxa_names(ps_genus) <- gsub("\\(|\\)", "", taxa_names(ps_genus))
#check
unique(tax_table(ps_genus)[,"Genus"] )## Taxonomy Table: [462 taxa by 1 taxonomic ranks]:
## Genus
## ASV_3181:g__VHS-B3-70 "VHS-B3-70"
## ASV_4916:g__Woesearchaeales "Woesearchaeales"
## ASV_1765:g__AKYG587 "AKYG587"
## ASV_1959:g__uncultured "uncultured"
## ASV_7171:g__Fimbriimonadaceae "Fimbriimonadaceae"
## ASV_2739:g__Opitutus "Opitutus"
## ASV_2126:g__Cloacibacterium "Cloacibacterium"
## ASV_4050:g__Subgroup_12 "Subgroup_12"
## ASV_3872:g__Hadarchaeales "Hadarchaeales"
## ASV_3878:g__Desulfuromonadaceae "Desulfuromonadaceae"
## ASV_1608:g__DEV008 "DEV008"
## ASV_5186:g__Sulfuritalea "Sulfuritalea"
## ASV_3275:g__Urania-1B-19_marine_sediment_group "Urania-1B-19_marine_sediment_group"
## ASV_840:g__WS2 "WS2"
## ASV_3640:g__Sediminispirochaeta "Sediminispirochaeta"
## ASV_3555:g__Anaerococcus "Anaerococcus"
## ASV_4206:g__Pelomonas "Pelomonas"
## ASV_98:g__Bacillus "Bacillus"
## ASV_702:g__Sulfuricurvum "Sulfuricurvum"
## ASV_1309:g__Comamonas "Comamonas"
## ASV_4053:g__Candidatus_Xiphinematobacter "Candidatus_Xiphinematobacter"
## ASV_5353:g__Amphiplicatus "Amphiplicatus"
## ASV_1850:g__Subgroup_23 "Subgroup_23"
## ASV_541:g__Thiobacillus "Thiobacillus"
## ASV_6812:g__Candidatus_Komeilibacteria "Candidatus_Komeilibacteria"
## ASV_6217:g__RB41 "RB41"
## ASV_2551:g__Curvibacter "Curvibacter"
## ASV_1166:g__Dadabacteriales "Dadabacteriales"
## ASV_511:g__Burkholderia-Caballeronia-Paraburkholderia "Burkholderia-Caballeronia-Paraburkholderia"
## ASV_2769:g__Kazania "Kazania"
## ASV_3171:g__UBA9983 "UBA9983"
## ASV_1838:g__Sphingomonas "Sphingomonas"
## ASV_2758:g__AKYH767 "AKYH767"
## ASV_5014:g__Pseudoclavibacter "Pseudoclavibacter"
## ASV_4921:g__UTBCD1 "UTBCD1"
## ASV_4553:g__[Cytophaga]_xylanolytica_group "[Cytophaga]_xylanolytica_group"
## ASV_3865:g__Kapabacteriales "Kapabacteriales"
## ASV_1920:g__CK-2C2-2 "CK-2C2-2"
## ASV_199:g__KF-JG30-B3 "KF-JG30-B3"
## ASV_1164:g__Chitinibacter "Chitinibacter"
## ASV_183:g__PLTA13 "PLTA13"
## ASV_5281:g__Dyadobacter "Dyadobacter"
## ASV_5:g__Aeromonas "Aeromonas"
## ASV_1384:g__DS-100 "DS-100"
## ASV_89:g__GIF9 "GIF9"
## ASV_3631:g__SCGC_AAA011-D5 "SCGC_AAA011-D5"
## ASV_5107:g__Succinivibrio "Succinivibrio"
## ASV_5917:g__Bacteroidales "Bacteroidales"
## ASV_5408:g__Dinghuibacter "Dinghuibacter"
## ASV_1880:g__Dechlorobacter "Dechlorobacter"
## ASV_1772:g__MB-A2-108 "MB-A2-108"
## ASV_1243:g__Acidibacter "Acidibacter"
## ASV_6811:g__Brevibacillus "Brevibacillus"
## ASV_730:g__LCP-89 "LCP-89"
## ASV_4530:g__Candidatus_Levybacteria "Candidatus_Levybacteria"
## ASV_4582:g__hgcI_clade "hgcI_clade"
## ASV_3225:g__Sphingobium "Sphingobium"
## ASV_3205:g__RBG-16-58-14 "RBG-16-58-14"
## ASV_146:g__Leptospirillum "Leptospirillum"
## ASV_5819:g__A4b "A4b"
## ASV_501:g__Proteus "Proteus"
## ASV_1065:g__SCGC_AB-179-E04 "SCGC_AB-179-E04"
## ASV_1288:g__Rhodoferax "Rhodoferax"
## ASV_1460:g__Anaeromyxobacter "Anaeromyxobacter"
## ASV_260:g__Subgroup_9 "Subgroup_9"
## ASV_1313:g__ADurb.Bin180 "ADurb.Bin180"
## ASV_997:g__Vermiphilaceae "Vermiphilaceae"
## ASV_6192:g__Lineage_IV "Lineage_IV"
## ASV_1909:g__KD4-96 "KD4-96"
## ASV_401:g__Calditrichaceae "Calditrichaceae"
## ASV_113:g__Zixibacteria "Zixibacteria"
## ASV_670:g__SC-I-84 "SC-I-84"
## ASV_495:g__RBG-13-46-9 "RBG-13-46-9"
## ASV_2598:g__10bav-F6 "10bav-F6"
## ASV_2545:g__37-13 "37-13"
## ASV_201:g__Entomoplasma "Entomoplasma"
## ASV_235:g__ADurb.Bin118 "ADurb.Bin118"
## ASV_470:g__Weissella "Weissella"
## ASV_393:g__TSAC18 "TSAC18"
## ASV_2272:g__S0134_terrestrial_group "S0134_terrestrial_group"
## ASV_84:g__Gemmobacter "Gemmobacter"
## ASV_502:g__Nitrospira "Nitrospira"
## ASV_139:g__661239 "661239"
## ASV_5789:g__Ochrobactrum "Ochrobactrum"
## ASV_383:g__Asteroleplasma "Asteroleplasma"
## ASV_2215:g__Pantoea "Pantoea"
## ASV_1264:g__Leptotrichia "Leptotrichia"
## ASV_4:g__Bifidobacterium "Bifidobacterium"
## ASV_5453:g__Chthoniobacter "Chthoniobacter"
## ASV_2401:g__Fictibacillus "Fictibacillus"
## ASV_2581:g__WOR-1 "WOR-1"
## ASV_122:g__Rheinheimera "Rheinheimera"
## ASV_2339:g__Aquabacterium "Aquabacterium"
## ASV_266:g__Hirschia "Hirschia"
## ASV_227:g__Staphylococcus "Staphylococcus"
## ASV_171:g__MSBL5 "MSBL5"
## ASV_5898:g__Agathobacter "Agathobacter"
## ASV_446:g__Clostridium_sensu_stricto_8 "Clostridium_sensu_stricto_8"
## ASV_1757:g__OLB14 "OLB14"
## ASV_4226:g__PS-B29 "PS-B29"
## ASV_97:g__Morganella "Morganella"
## ASV_5731:g__Obscuribacteraceae "Obscuribacteraceae"
## ASV_296:g__ER-E4-17 "ER-E4-17"
## ASV_2238:g__Chryseolinea "Chryseolinea"
## ASV_587:g__Subgroup_2 "Subgroup_2"
## ASV_5651:g__Sulfuriferula "Sulfuriferula"
## ASV_886:g__Dechloromonas "Dechloromonas"
## ASV_2301:g__Streptomyces "Streptomyces"
## ASV_1093:g__Providencia "Providencia"
## ASV_1377:g__Candidatus_Collierbacteria "Candidatus_Collierbacteria"
## ASV_1646:g__UCG-004 "UCG-004"
## ASV_596:g__Subgroup_17 "Subgroup_17"
## ASV_352:g__mle1-7 "mle1-7"
## ASV_134:g__Desulfobacca "Desulfobacca"
## ASV_1992:g__AKAU3564_sediment_group "AKAU3564_sediment_group"
## ASV_1091:g__MBMPE27 "MBMPE27"
## ASV_1785:g__Pajaroellobacter "Pajaroellobacter"
## ASV_463:g__SJA-28 "SJA-28"
## ASV_111:g__MBNT15 "MBNT15"
## ASV_276:g__BSV26 "BSV26"
## ASV_6374:g__Desulfatirhabdium "Desulfatirhabdium"
## ASV_59:g__Pectinatus "Pectinatus"
## ASV_563:g__GOUTB8 "GOUTB8"
## ASV_5086:g__Acidobacteriota "Acidobacteriota"
## ASV_4953:g__Thermodesulfovibrio "Thermodesulfovibrio"
## ASV_1516:g__Kineosporia "Kineosporia"
## ASV_58:g__Bradyrhizobium "Bradyrhizobium"
## ASV_683:g__Syntrophorhabdus "Syntrophorhabdus"
## ASV_1030:g__bacteriap25 "bacteriap25"
## ASV_414:g__SAR202_clade "SAR202_clade"
## ASV_7033:g__Micromonospora "Micromonospora"
## ASV_2640:g__SJA-15 "SJA-15"
## ASV_3609:g__Paucibacter "Paucibacter"
## ASV_3740:g__Brevundimonas "Brevundimonas"
## ASV_3379:g__Lawsonella "Lawsonella"
## ASV_2528:g__Marinimicrobia_SAR406_clade "Marinimicrobia_(SAR406_clade)"
## ASV_516:g__Nocardiopsis "Nocardiopsis"
## ASV_4827:g__Candidatus_Doudnabacteria "Candidatus_Doudnabacteria"
## ASV_126:g__H3.93 "H3.93"
## ASV_5109:g__Gemmata "Gemmata"
## ASV_1401:g__PAUC26f "PAUC26f"
## ASV_1936:g__Cupriavidus "Cupriavidus"
## ASV_4746:g__Treponema "Treponema"
## ASV_644:g__GAL15 "GAL15"
## ASV_1967:g__Candidatus_Kaiserbacteria "Candidatus_Kaiserbacteria"
## ASV_4485:g__Vagococcus "Vagococcus"
## ASV_204:g__Methylobacterium-Methylorubrum "Methylobacterium-Methylorubrum"
## ASV_3596:g__Candidatus_Moranbacteria "Candidatus_Moranbacteria"
## ASV_1903:g__Z195MB87 "Z195MB87"
## ASV_350:g__Clostridium_sensu_stricto_5 "Clostridium_sensu_stricto_5"
## ASV_3246:g__IMCC26256 "IMCC26256"
## ASV_877:g__Flavobacterium "Flavobacterium"
## ASV_67:g__Spirochaeta "Spirochaeta"
## ASV_92:g__GIF3 "GIF3"
## ASV_465:g__Sh765B-TzT-20 "Sh765B-TzT-20"
## ASV_2317:g__Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium "Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium"
## ASV_6839:g__P3OB-42 "P3OB-42"
## ASV_599:g__Chitinimonas "Chitinimonas"
## ASV_1092:g__Marine_Benthic_Group_D_and_DHVEG-1 "Marine_Benthic_Group_D_and_DHVEG-1"
## ASV_1392:g__Romboutsia "Romboutsia"
## ASV_1256:g__Sericytochromatia "Sericytochromatia"
## ASV_3182:g__Ga0077536 "Ga0077536"
## ASV_148:g__S085 "S085"
## ASV_226:g__Streptococcus "Streptococcus"
## ASV_429:g__Hypnocyclicus "Hypnocyclicus"
## ASV_356:g__GOUTA6 "GOUTA6"
## ASV_6978:g__Hafnia-Obesumbacterium "Hafnia-Obesumbacterium"
## ASV_6342:g__Gracilibacteria "Gracilibacteria"
## ASV_32:g__Aminicenantales "Aminicenantales"
## ASV_1238:g__SPG12-343-353-B75 "SPG12-343-353-B75"
## ASV_2547:g__Lacunisphaera "Lacunisphaera"
## ASV_769:g__PHOS-HE36 "PHOS-HE36"
## ASV_4642:g__MSB-5E12 "MSB-5E12"
## ASV_5286:g__Dyella "Dyella"
## ASV_362:g__Dysgonomonas "Dysgonomonas"
## ASV_1519:g__Subgroup_22 "Subgroup_22"
## ASV_132:g__SB-5 "SB-5"
## ASV_6145:g__Pedosphaera "Pedosphaera"
## ASV_1140:g__AKIW659 "AKIW659"
## ASV_925:g__SEEP-SRB1 "SEEP-SRB1"
## ASV_2297:g__Schlegelella "Schlegelella"
## ASV_5112:g__Fimbriiglobus "Fimbriiglobus"
## ASV_6393:g__Phyllobacterium "Phyllobacterium"
## ASV_852:g__FS117-23B-02 "FS117-23B-02"
## ASV_3317:g__Gitt-GS-136 "Gitt-GS-136"
## ASV_2927:g__Hydrogenedensaceae "Hydrogenedensaceae"
## ASV_5128:g__Dehalogenimonas "Dehalogenimonas"
## ASV_1637:g__LD1-PA32 "LD1-PA32"
## ASV_1492:g__CCM11a "CCM11a"
## ASV_1220:g__CG2-30-50-142 "CG2-30-50-142"
## ASV_6132:g__Fusobacterium "Fusobacterium"
## ASV_160:g__KF-JG30-C25 "KF-JG30-C25"
## ASV_1308:g__S-BQ2-57_soil_group "S-BQ2-57_soil_group"
## ASV_3908:g__LWQ8 "LWQ8"
## ASV_223:g__BD2-11_terrestrial_group "BD2-11_terrestrial_group"
## ASV_2414:g__Deferrisoma "Deferrisoma"
## ASV_21:g__Micrococcus "Micrococcus"
## ASV_4595:g__Clostridium_sensu_stricto_9 "Clostridium_sensu_stricto_9"
## ASV_585:g__Subgroup_7 "Subgroup_7"
## ASV_1851:g__BSN166 "BSN166"
## ASV_1702:g__WD2101_soil_group "WD2101_soil_group"
## ASV_297:g__Calorithrix "Calorithrix"
## ASV_5723:g__Pir1_lineage "Pir1_lineage"
## ASV_4999:g__Run-SP154 "Run-SP154"
## ASV_1564:g__Thioalkalispira-Sulfurivermis "Thioalkalispira-Sulfurivermis"
## ASV_4344:g__OC31 "OC31"
## ASV_317:g__Megasphaera "Megasphaera"
## ASV_6814:g__CI75cm.2.12 "CI75cm.2.12"
## ASV_27:g__TRA3-20 "TRA3-20"
## ASV_921:g__Oceanobacillus "Oceanobacillus"
## ASV_1124:g__Faecalibacterium "Faecalibacterium"
## ASV_1990:g__Candidatus_Adlerbacteria "Candidatus_Adlerbacteria"
## ASV_834:g__Thermoanaerobaculum "Thermoanaerobaculum"
## ASV_5403:g__Sphingobacterium "Sphingobacterium"
## ASV_320:g__EPR3968-O8a-Bc78 "EPR3968-O8a-Bc78"
## ASV_1009:g__Moduliflexaceae "Moduliflexaceae"
## ASV_1932:g__Candidatus_Peregrinibacteria "Candidatus_Peregrinibacteria"
## ASV_1075:g__Rhodococcus "Rhodococcus"
## ASV_1491:g__mle1-8 "mle1-8"
## ASV_167:g__Deinococcus "Deinococcus"
## ASV_4158:g__WN-HWB-116 "WN-HWB-116"
## ASV_1873:g__Stenotrophomonas "Stenotrophomonas"
## ASV_4902:g__JdFR-54 "JdFR-54"
## ASV_802:g__Raoultella "Raoultella"
## ASV_1249:g__Crenobacter "Crenobacter"
## ASV_1520:g__WX65 "WX65"
## ASV_131:g__A0839 "A0839"
## ASV_221:g__vadinBA26 "vadinBA26"
## ASV_1950:g__Acetothermiia "Acetothermiia"
## ASV_1518:g__Shimwellia "Shimwellia"
## ASV_693:g__RCP2-54 "RCP2-54"
## ASV_4558:g__Luteolibacter "Luteolibacter"
## ASV_5802:g__Desulfovibrio "Desulfovibrio"
## ASV_2653:g__Pediococcus "Pediococcus"
## ASV_325:g__WWE3 "WWE3"
## ASV_3582:g__Schekmanbacteria "Schekmanbacteria"
## ASV_1341:g__Maricaulis "Maricaulis"
## ASV_2727:g__ADurb.Bin063-1 "ADurb.Bin063-1"
## ASV_1779:g__Caldithrix "Caldithrix"
## ASV_194:g__FW22 "FW22"
## ASV_1254:g__Arcobacter "Arcobacter"
## ASV_1684:g__Fimbriimonadales "Fimbriimonadales"
## ASV_6678:g__Lacibacter "Lacibacter"
## ASV_3700:g__B103G10 "B103G10"
## ASV_1919:g__Pir4_lineage "Pir4_lineage"
## ASV_2957:g__Pedosphaeraceae "Pedosphaeraceae"
## ASV_5692:g__Nordella "Nordella"
## ASV_286:g__Lactococcus "Lactococcus"
## ASV_5697:g__Steroidobacter "Steroidobacter"
## ASV_3155:g__Z-35 "Z-35"
## ASV_4374:g__Moheibacter "Moheibacter"
## ASV_1691:g__Desulfomonile "Desulfomonile"
## ASV_2342:g__Syntrophobacter "Syntrophobacter"
## ASV_26:g__Napoli-4B-65 "Napoli-4B-65"
## ASV_2928:g__Sandaracinus "Sandaracinus"
## ASV_1818:g__Domibacillus "Domibacillus"
## ASV_705:g__Bauldia "Bauldia"
## ASV_24:g__Escherichia-Shigella "Escherichia-Shigella"
## ASV_2434:g__Pla4_lineage "Pla4_lineage"
## ASV_5301:g__Megamonas "Megamonas"
## ASV_4731:g__Actinomyces "Actinomyces"
## ASV_118:g__MND1 "MND1"
## ASV_1862:g__Subgroup_25 "Subgroup_25"
## ASV_2521:g__Candidatus_Liptonbacteria "Candidatus_Liptonbacteria"
## ASV_6258:g__WS4 "WS4"
## ASV_477:g__Subgroup_15 "Subgroup_15"
## ASV_1615:g__B2M28 "B2M28"
## ASV_2747:g__Cutibacterium "Cutibacterium"
## ASV_6917:g__Chitinophaga "Chitinophaga"
## ASV_4295:g__Subgroup_10 "Subgroup_10"
## ASV_3212:g__BIrii41 "BIrii41"
## ASV_1276:g__Gammaproteobacteria "Gammaproteobacteria"
## ASV_208:g__Dehalococcoidia "Dehalococcoidia"
## ASV_5788:g__Candidatus_Uhrbacteria "Candidatus_Uhrbacteria"
## ASV_7160:g__Pla3_lineage "Pla3_lineage"
## ASV_3495:g__Blfdi19 "Blfdi19"
## ASV_1489:g__Candidatus_Magasanikbacteria "Candidatus_Magasanikbacteria"
## ASV_5104:g__Thiogranum "Thiogranum"
## ASV_1027:g__ODP1230B23.02 "ODP1230B23.02"
## ASV_4874:g__Eel-36e1D6 "Eel-36e1D6"
## ASV_115:g__Bathyarchaeia "Bathyarchaeia"
## ASV_1325:g__Terrimonas "Terrimonas"
## ASV_2644:g__Clostridium_sensu_stricto_13 "Clostridium_sensu_stricto_13"
## ASV_105:g__Sva0081_sediment_group "Sva0081_sediment_group"
## ASV_7:g__Sva0485 "Sva0485"
## ASV_137:g__Exiguobacterium "Exiguobacterium"
## ASV_374:g__Thermoflexus "Thermoflexus"
## ASV_3438:g__oc32 "oc32"
## ASV_1517:g__JG36-GS-52 "JG36-GS-52"
## ASV_3269:g__Gemmatimonadaceae "Gemmatimonadaceae"
## ASV_678:g__RBG-16-49-21 "RBG-16-49-21"
## ASV_5555:g__Shewanella "Shewanella"
## ASV_6776:g__Bradymonadales "Bradymonadales"
## ASV_22:g__Lactobacillus "Lactobacillus"
## ASV_3261:g__SH-PL14 "SH-PL14"
## ASV_192:g__4-29-1 "4-29-1"
## ASV_709:g__Subgroup_13 "Subgroup_13"
## ASV_3165:g__Rhodomicrobium "Rhodomicrobium"
## ASV_738:g__Bacteroidetes_BD2-2 "Bacteroidetes_BD2-2"
## ASV_456:g__SBR1031 "SBR1031"
## ASV_1642:g__Candidatus_Paceibacter "Candidatus_Paceibacter"
## ASV_6714:g__Proteiniphilum "Proteiniphilum"
## ASV_593:g__Sporacetigenium "Sporacetigenium"
## ASV_198:g__Serratia "Serratia"
## ASV_1797:g__Odinarchaeia "Odinarchaeia"
## ASV_4836:g__Candidatus_Alysiosphaera "Candidatus_Alysiosphaera"
## ASV_3711:g__Chitinilyticum "Chitinilyticum"
## ASV_6655:g__Subdoligranulum "Subdoligranulum"
## ASV_3402:g__Ruminiclostridium "Ruminiclostridium"
## ASV_3089:g__Altererythrobacter "Altererythrobacter"
## ASV_4278:g__LD-RB-34 "LD-RB-34"
## ASV_200:g__MD2902-B12 "MD2902-B12"
## ASV_6053:g__Holophagaceae "Holophagaceae"
## ASV_6512:g__UCG-012 "UCG-012"
## ASV_5879:g__Ferruginibacter "Ferruginibacter"
## ASV_2070:g__SG8-4 "SG8-4"
## ASV_55:g__Acinetobacter "Acinetobacter"
## ASV_272:g__Pedomicrobium "Pedomicrobium"
## ASV_51:g__Enhydrobacter "Enhydrobacter"
## ASV_459:g__Subgroup_21 "Subgroup_21"
## ASV_3096:g__KI89A_clade "KI89A_clade"
## ASV_5622:g__Fervidicella "Fervidicella"
## ASV_4617:g__Candidatus_Chisholmbacteria "Candidatus_Chisholmbacteria"
## ASV_265:g__Latescibacteraceae "Latescibacteraceae"
## ASV_1011:g__DG-56 "DG-56"
## ASV_110:g__ZOR0006 "ZOR0006"
## ASV_4811:g__Candidatus_Thiosymbion "Candidatus_Thiosymbion"
## ASV_3791:g__A21b "A21b"
## ASV_1189:g__Luteimonas "Luteimonas"
## ASV_240:g__Sh765B-AG-111 "Sh765B-AG-111"
## ASV_990:g__Dongia "Dongia"
## ASV_6472:g__Rivicola "Rivicola"
## ASV_2371:g__SHA-26 "SHA-26"
## ASV_2534:g__Tumebacillus "Tumebacillus"
## ASV_1099:g__Hyphomicrobium "Hyphomicrobium"
## ASV_2232:g__Inhella "Inhella"
## ASV_3399:g__P9X2b3D02 "P9X2b3D02"
## ASV_2599:g__Candidatus_Latescibacter "Candidatus_Latescibacter"
## ASV_2105:g__Candidatus_Solibacter "Candidatus_Solibacter"
## ASV_4958:g__Woeseia "Woeseia"
## ASV_5875:g__Entotheonellaceae "Entotheonellaceae"
## ASV_1465:g__Clostridium_sensu_stricto_6 "Clostridium_sensu_stricto_6"
## ASV_1991:g__Subgroup_5 "Subgroup_5"
## ASV_4177:g__Limibaculum "Limibaculum"
## ASV_5483:g__[Polyangium]_brachysporum_group "[Polyangium]_brachysporum_group"
## ASV_73:g__t0.6.f "t0.6.f"
## ASV_69:g__Desulfatiglans "Desulfatiglans"
## ASV_2235:g__Ensifer "Ensifer"
## ASV_5272:g__Candidatus_Portnoybacteria "Candidatus_Portnoybacteria"
## ASV_5890:g__Candidatus_Methylospira "Candidatus_Methylospira"
## ASV_4434:g__Nitrolancea "Nitrolancea"
## ASV_4377:g__Pla1_lineage "Pla1_lineage"
## ASV_360:g__NB1-j "NB1-j"
## ASV_70:g__SCGC-AB-539-J10 "SCGC-AB-539-J10"
## ASV_6034:g__SAR324_cladeMarine_group_B "SAR324_clade(Marine_group_B)"
## ASV_1291:g__Parcubacteria "Parcubacteria"
## ASV_1461:g__vadinHA49 "vadinHA49"
## ASV_2782:g__Saccharimonadales "Saccharimonadales"
## ASV_433:g__Latescibacterota "Latescibacterota"
## ASV_1706:g__AKAU4049 "AKAU4049"
## ASV_4985:g__Novosphingobium "Novosphingobium"
## ASV_5496:g__EC3 "EC3"
## ASV_3068:g__WCHB1-81 "WCHB1-81"
## ASV_141:g__Corynebacterium "Corynebacterium"
## ASV_676:g__Clostridium_sensu_stricto_1 "Clostridium_sensu_stricto_1"
## ASV_2785:g__Haliangium "Haliangium"
## ASV_4788:g__Pirellula "Pirellula"
## ASV_4272:g__Puia "Puia"
## ASV_3:g__Enterobacter "Enterobacter"
## ASV_480:g__TM7a "TM7a"
## ASV_5877:g__Solitalea "Solitalea"
## ASV_313:g__Ellin516 "Ellin516"
## ASV_3257:g__Candidatus_Woesebacteria "Candidatus_Woesebacteria"
## ASV_39:g__Sulfurifustis "Sulfurifustis"
## ASV_450:g__Brachybacterium "Brachybacterium"
## ASV_5386:g__Permianibacter "Permianibacter"
## ASV_4805:g__3PJM14 "3PJM14"
## ASV_5535:g__Leeia "Leeia"
## ASV_519:g__Poribacteria "Poribacteria"
## ASV_3723:g__Amsterdam-1B-07 "Amsterdam-1B-07"
## ASV_34:g__Subgroup_18 "Subgroup_18"
## ASV_2077:g__Lysobacter "Lysobacter"
## ASV_3408:g__wb1-A12 "wb1-A12"
## ASV_228:g__Ignavibacterium "Ignavibacterium"
## ASV_539:g__Pelolinea "Pelolinea"
## ASV_4701:g__UASB-TL25 "UASB-TL25"
## ASV_5499:g__S-70 "S-70"
## ASV_7088:g__EF100-94H03 "EF100-94H03"
## ASV_6730:g__Candidatus_Nomurabacteria "Candidatus_Nomurabacteria"
## ASV_5459:g__MM2 "MM2"
## ASV_1287:g__BD1-7_clade "BD1-7_clade"
## ASV_1:g__Enterococcus "Enterococcus"
## ASV_4188:g__WPS-2 "WPS-2"
## ASV_4223:g__Parapedobacter "Parapedobacter"
## ASV_112:g__Rokubacteriales "Rokubacteriales"
## ASV_857:g__Candidatus_Yanofskybacteria "Candidatus_Yanofskybacteria"
## ASV_5663:g__966-1 "966-1"
## ASV_775:g__Bacteroidetes_vadinHA17 "Bacteroidetes_vadinHA17"
## ASV_1356:g__Ralstonia "Ralstonia"
## ASV_819:g__RBG-13-54-9 "RBG-13-54-9"
## ASV_3952:g__Subgroup_11 "Subgroup_11"
## ASV_746:g__Massilia "Massilia"
## ASV_4792:g__Verrucomicrobium "Verrucomicrobium"
## ASV_6367:g__Candidatus_Jorgensenbacteria "Candidatus_Jorgensenbacteria"
## ASV_3672:g__OM190 "OM190"
## ASV_130:g__Acetobacter "Acetobacter"
## ASV_4124:g__Sideroxydans "Sideroxydans"
## ASV_4111:g__Tolumonas "Tolumonas"
## ASV_752:g__Bryobacter "Bryobacter"
## ASV_703:g__Candidatus_Omnitrophus "Candidatus_Omnitrophus"
## ASV_3320:g__Rhodobacter "Rhodobacter"
## ASV_5092:g__Longimicrobiaceae "Longimicrobiaceae"
## ASV_3517:g__OPB41 "OPB41"
## ASV_66:g__Sh765B-TzT-35 "Sh765B-TzT-35"
## ASV_4666:g__CL500-29_marine_group "CL500-29_marine_group"
## ASV_72:g__Vogesella "Vogesella"
## ASV_3529:g__Thermincola "Thermincola"
## ASV_608:g__CCD24 "CCD24"
## ASV_590:g__Thiobacter "Thiobacter"
## ASV_5373:g__Microgenomatia "Microgenomatia"
## ASV_5790:g__Bacteroidetes_VC2.1_Bac22 "Bacteroidetes_VC2.1_Bac22"
## ASV_41:g__B1-7BS "B1-7BS"
## ASV_942:g__Olsenella "Olsenella"
## ASV_6267:g__Devosia "Devosia"
## ASV_725:g__JG30-KF-CM66 "JG30-KF-CM66"
## ASV_513:g__Sumerlaea "Sumerlaea"
## ASV_2511:g__Reyranella "Reyranella"
## ASV_4756:g__Blastococcus "Blastococcus"
## ASV_1041:g__Mesorhizobium "Mesorhizobium"
## ASV_5164:g__RBG-16-55-12 "RBG-16-55-12"
## ASV_1001:g__AB-539-J10 "AB-539-J10"
## ASV_555:g__Hydrogenophaga "Hydrogenophaga"
## ASV_856:g__Pseudolabrys "Pseudolabrys"
## ASV_184:g__Vicinamibacteraceae "Vicinamibacteraceae"
## ASV_2719:g__SM23-31 "SM23-31"
## ASV_4463:g__OM27_clade "OM27_clade"
## ASV_797:g__Nitratiruptor "Nitratiruptor"
## ASV_1078:g__GWE2-31-10 "GWE2-31-10"
## ASV_3338:g__Paenarthrobacter "Paenarthrobacter"
## ASV_6583:g__DEV114 "DEV114"
## ASV_8:g__Prevotella "Prevotella"
## ASV_3311:g__Nitrosopumilaceae "Nitrosopumilaceae"
## ASV_2322:g__Methyloligellaceae "Methyloligellaceae"
## ASV_1839:g__Leuconostoc "Leuconostoc"
## ASV_5937:g__Azohydromonas "Azohydromonas"
## ASV_1423:g__Gaiella "Gaiella"
## ASV_2282:g__Vulgatibacter "Vulgatibacter"
## ASV_1114:g__Ellin6067 "Ellin6067"
## ASV_1503:g__MSB-5B2 "MSB-5B2"
## ASV_3486:g__Lentimicrobiaceae "Lentimicrobiaceae"
## ASV_812:g__Z114MB74 "Z114MB74"
## ASV_979:g__Edwardsbacteria "Edwardsbacteria"
## ASV_611:g__11-24 "11-24"
## ASV_7101:g__JG36-TzT-191 "JG36-TzT-191"
## ASV_6:g__Pseudomonas "Pseudomonas"
## ASV_935:g__Dialister "Dialister"
## ASV_2933:g__Pseudochrobactrum "Pseudochrobactrum"
## ASV_328:g__TA06 "TA06"
## ASV_355:g__SWB02 "SWB02"
## ASV_1639:g__Fluviicola "Fluviicola"
## ASV_2616:g__Kytococcus "Kytococcus"
## ASV_1246:g__Delftia "Delftia"
# check data
table(meta(ps)$Site, meta(ps)$Depth_cm)##
## 0-15 16-30 31-45 50-75
## Site_1 2 2 2 2
## Site_2 2 2 2 2
## Site_3 2 2 2 2
## Site_4 1 1 1 1
## Site_5 2 2 2 2
## Site_6 2 2 2 0
## Site_7 2 2 2 2
#subsets
ps_s1 <- subset_samples(ps, Site == "Site_1")
ps_s2 <- subset_samples(ps, Site == "Site_2")
ps_s3 <- subset_samples(ps, Site == "Site_3")
ps_s4 <- subset_samples(ps, Site == "Site_4")
ps_s5 <- subset_samples(ps, Site == "Site_5")
ps_s6 <- subset_samples(ps, Site == "Site_6")
ps_s7 <- subset_samples(ps, Site == "Site_7")
#show
ps_s1## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 7097 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 7097 taxa by 7 taxonomic ranks ]
ps_s2## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 7097 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 7097 taxa by 7 taxonomic ranks ]
ps_s3## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 7097 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 7097 taxa by 7 taxonomic ranks ]
ps_s4## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 7097 taxa and 4 samples ]
## sample_data() Sample Data: [ 4 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 7097 taxa by 7 taxonomic ranks ]
ps_s5## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 7097 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 7097 taxa by 7 taxonomic ranks ]
ps_s6## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 7097 taxa and 6 samples ]
## sample_data() Sample Data: [ 6 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 7097 taxa by 7 taxonomic ranks ]
ps_s7## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 7097 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 7097 taxa by 7 taxonomic ranks ]
#check data
microbiome::summarize_phyloseq(ps_s1)## Compositional = NO2
## 1] Min. number of reads = 1486352] Max. number of reads = 2188773] Total number of reads = 14291624] Average number of reads = 178645.255] Median number of reads = 1712277] Sparsity = 0.5537374947160776] Any OTU sum to 1 or less? YES8] Number of singletons = 7879] Percent of OTUs that are singletons
## (i.e. exactly one read detected across all samples)010] Number of sample variables are: 29SampleIDDepthDepth_cmSiteSample_typeDateLocationTemp...C..intSalnidad..ppm.intSalinidad..ppt.intpH.intpHsupRedox.medido.en.campo1Redox.mV..Eh.intTemperatura...C.supSalnidad..ppm.supSalinidad..ppt.supRedox.medido.en.campoRedox.mV..Eh.supS.2..mg.l..campoml.de.muestra1Factor.de.dilucion1S.2..mg.l.ID_MSO.....mg.l..campoml.de.muestra2Factor.de.dilucion2SO.....mg.l.2
## [[1]]
## [1] "1] Min. number of reads = 148635"
##
## [[2]]
## [1] "2] Max. number of reads = 218877"
##
## [[3]]
## [1] "3] Total number of reads = 1429162"
##
## [[4]]
## [1] "4] Average number of reads = 178645.25"
##
## [[5]]
## [1] "5] Median number of reads = 171227"
##
## [[6]]
## [1] "7] Sparsity = 0.553737494716077"
##
## [[7]]
## [1] "6] Any OTU sum to 1 or less? YES"
##
## [[8]]
## [1] "8] Number of singletons = 787"
##
## [[9]]
## [1] "9] Percent of OTUs that are singletons \n (i.e. exactly one read detected across all samples)0"
##
## [[10]]
## [1] "10] Number of sample variables are: 29"
##
## [[11]]
## [1] "Sample" "ID" "Depth"
## [4] "Depth_cm" "Site" "Sample_type"
## [7] "Date" "Location" "Temp...C..int"
## [10] "Salnidad..ppm.int" "Salinidad..ppt.int" "pH.int"
## [13] "pHsup" "Redox.medido.en.campo1" "Redox.mV..Eh.int"
## [16] "Temperatura...C.sup" "Salnidad..ppm.sup" "Salinidad..ppt.sup"
## [19] "Redox.medido.en.campo" "Redox.mV..Eh.sup" "S.2..mg.l..campo"
## [22] "ml.de.muestra1" "Factor.de.dilucion1" "S.2..mg.l."
## [25] "ID_M" "SO.....mg.l..campo" "ml.de.muestra2"
## [28] "Factor.de.dilucion2" "SO.....mg.l."
microbiome::summarize_phyloseq(ps_s2)## Compositional = NO2
## 1] Min. number of reads = 1323482] Max. number of reads = 2205473] Total number of reads = 13798474] Average number of reads = 172480.8755] Median number of reads = 173469.57] Sparsity = 0.5038572636325216] Any OTU sum to 1 or less? YES8] Number of singletons = 12669] Percent of OTUs that are singletons
## (i.e. exactly one read detected across all samples)010] Number of sample variables are: 29SampleIDDepthDepth_cmSiteSample_typeDateLocationTemp...C..intSalnidad..ppm.intSalinidad..ppt.intpH.intpHsupRedox.medido.en.campo1Redox.mV..Eh.intTemperatura...C.supSalnidad..ppm.supSalinidad..ppt.supRedox.medido.en.campoRedox.mV..Eh.supS.2..mg.l..campoml.de.muestra1Factor.de.dilucion1S.2..mg.l.ID_MSO.....mg.l..campoml.de.muestra2Factor.de.dilucion2SO.....mg.l.2
## [[1]]
## [1] "1] Min. number of reads = 132348"
##
## [[2]]
## [1] "2] Max. number of reads = 220547"
##
## [[3]]
## [1] "3] Total number of reads = 1379847"
##
## [[4]]
## [1] "4] Average number of reads = 172480.875"
##
## [[5]]
## [1] "5] Median number of reads = 173469.5"
##
## [[6]]
## [1] "7] Sparsity = 0.503857263632521"
##
## [[7]]
## [1] "6] Any OTU sum to 1 or less? YES"
##
## [[8]]
## [1] "8] Number of singletons = 1266"
##
## [[9]]
## [1] "9] Percent of OTUs that are singletons \n (i.e. exactly one read detected across all samples)0"
##
## [[10]]
## [1] "10] Number of sample variables are: 29"
##
## [[11]]
## [1] "Sample" "ID" "Depth"
## [4] "Depth_cm" "Site" "Sample_type"
## [7] "Date" "Location" "Temp...C..int"
## [10] "Salnidad..ppm.int" "Salinidad..ppt.int" "pH.int"
## [13] "pHsup" "Redox.medido.en.campo1" "Redox.mV..Eh.int"
## [16] "Temperatura...C.sup" "Salnidad..ppm.sup" "Salinidad..ppt.sup"
## [19] "Redox.medido.en.campo" "Redox.mV..Eh.sup" "S.2..mg.l..campo"
## [22] "ml.de.muestra1" "Factor.de.dilucion1" "S.2..mg.l."
## [25] "ID_M" "SO.....mg.l..campo" "ml.de.muestra2"
## [28] "Factor.de.dilucion2" "SO.....mg.l."
microbiome::summarize_phyloseq(ps_s3)## Compositional = NO2
## 1] Min. number of reads = 1019792] Max. number of reads = 6748463] Total number of reads = 16868954] Average number of reads = 210861.8755] Median number of reads = 1415347] Sparsity = 0.6287691982527836] Any OTU sum to 1 or less? YES8] Number of singletons = 23169] Percent of OTUs that are singletons
## (i.e. exactly one read detected across all samples)010] Number of sample variables are: 29SampleIDDepthDepth_cmSiteSample_typeDateLocationTemp...C..intSalnidad..ppm.intSalinidad..ppt.intpH.intpHsupRedox.medido.en.campo1Redox.mV..Eh.intTemperatura...C.supSalnidad..ppm.supSalinidad..ppt.supRedox.medido.en.campoRedox.mV..Eh.supS.2..mg.l..campoml.de.muestra1Factor.de.dilucion1S.2..mg.l.ID_MSO.....mg.l..campoml.de.muestra2Factor.de.dilucion2SO.....mg.l.2
## [[1]]
## [1] "1] Min. number of reads = 101979"
##
## [[2]]
## [1] "2] Max. number of reads = 674846"
##
## [[3]]
## [1] "3] Total number of reads = 1686895"
##
## [[4]]
## [1] "4] Average number of reads = 210861.875"
##
## [[5]]
## [1] "5] Median number of reads = 141534"
##
## [[6]]
## [1] "7] Sparsity = 0.628769198252783"
##
## [[7]]
## [1] "6] Any OTU sum to 1 or less? YES"
##
## [[8]]
## [1] "8] Number of singletons = 2316"
##
## [[9]]
## [1] "9] Percent of OTUs that are singletons \n (i.e. exactly one read detected across all samples)0"
##
## [[10]]
## [1] "10] Number of sample variables are: 29"
##
## [[11]]
## [1] "Sample" "ID" "Depth"
## [4] "Depth_cm" "Site" "Sample_type"
## [7] "Date" "Location" "Temp...C..int"
## [10] "Salnidad..ppm.int" "Salinidad..ppt.int" "pH.int"
## [13] "pHsup" "Redox.medido.en.campo1" "Redox.mV..Eh.int"
## [16] "Temperatura...C.sup" "Salnidad..ppm.sup" "Salinidad..ppt.sup"
## [19] "Redox.medido.en.campo" "Redox.mV..Eh.sup" "S.2..mg.l..campo"
## [22] "ml.de.muestra1" "Factor.de.dilucion1" "S.2..mg.l."
## [25] "ID_M" "SO.....mg.l..campo" "ml.de.muestra2"
## [28] "Factor.de.dilucion2" "SO.....mg.l."
microbiome::summarize_phyloseq(ps_s4)## Compositional = NO2
## 1] Min. number of reads = 1276312] Max. number of reads = 2065043] Total number of reads = 7156094] Average number of reads = 178902.255] Median number of reads = 1907377] Sparsity = 0.4716077215724956] Any OTU sum to 1 or less? YES8] Number of singletons = 16899] Percent of OTUs that are singletons
## (i.e. exactly one read detected across all samples)010] Number of sample variables are: 29SampleIDDepthDepth_cmSiteSample_typeDateLocationTemp...C..intSalnidad..ppm.intSalinidad..ppt.intpH.intpHsupRedox.medido.en.campo1Redox.mV..Eh.intTemperatura...C.supSalnidad..ppm.supSalinidad..ppt.supRedox.medido.en.campoRedox.mV..Eh.supS.2..mg.l..campoml.de.muestra1Factor.de.dilucion1S.2..mg.l.ID_MSO.....mg.l..campoml.de.muestra2Factor.de.dilucion2SO.....mg.l.2
## [[1]]
## [1] "1] Min. number of reads = 127631"
##
## [[2]]
## [1] "2] Max. number of reads = 206504"
##
## [[3]]
## [1] "3] Total number of reads = 715609"
##
## [[4]]
## [1] "4] Average number of reads = 178902.25"
##
## [[5]]
## [1] "5] Median number of reads = 190737"
##
## [[6]]
## [1] "7] Sparsity = 0.471607721572495"
##
## [[7]]
## [1] "6] Any OTU sum to 1 or less? YES"
##
## [[8]]
## [1] "8] Number of singletons = 1689"
##
## [[9]]
## [1] "9] Percent of OTUs that are singletons \n (i.e. exactly one read detected across all samples)0"
##
## [[10]]
## [1] "10] Number of sample variables are: 29"
##
## [[11]]
## [1] "Sample" "ID" "Depth"
## [4] "Depth_cm" "Site" "Sample_type"
## [7] "Date" "Location" "Temp...C..int"
## [10] "Salnidad..ppm.int" "Salinidad..ppt.int" "pH.int"
## [13] "pHsup" "Redox.medido.en.campo1" "Redox.mV..Eh.int"
## [16] "Temperatura...C.sup" "Salnidad..ppm.sup" "Salinidad..ppt.sup"
## [19] "Redox.medido.en.campo" "Redox.mV..Eh.sup" "S.2..mg.l..campo"
## [22] "ml.de.muestra1" "Factor.de.dilucion1" "S.2..mg.l."
## [25] "ID_M" "SO.....mg.l..campo" "ml.de.muestra2"
## [28] "Factor.de.dilucion2" "SO.....mg.l."
microbiome::summarize_phyloseq(ps_s5)## Compositional = NO2
## 1] Min. number of reads = 970812] Max. number of reads = 2982253] Total number of reads = 14629924] Average number of reads = 1828745] Median number of reads = 1679587] Sparsity = 0.5306643652247436] Any OTU sum to 1 or less? YES8] Number of singletons = 8929] Percent of OTUs that are singletons
## (i.e. exactly one read detected across all samples)010] Number of sample variables are: 29SampleIDDepthDepth_cmSiteSample_typeDateLocationTemp...C..intSalnidad..ppm.intSalinidad..ppt.intpH.intpHsupRedox.medido.en.campo1Redox.mV..Eh.intTemperatura...C.supSalnidad..ppm.supSalinidad..ppt.supRedox.medido.en.campoRedox.mV..Eh.supS.2..mg.l..campoml.de.muestra1Factor.de.dilucion1S.2..mg.l.ID_MSO.....mg.l..campoml.de.muestra2Factor.de.dilucion2SO.....mg.l.2
## [[1]]
## [1] "1] Min. number of reads = 97081"
##
## [[2]]
## [1] "2] Max. number of reads = 298225"
##
## [[3]]
## [1] "3] Total number of reads = 1462992"
##
## [[4]]
## [1] "4] Average number of reads = 182874"
##
## [[5]]
## [1] "5] Median number of reads = 167958"
##
## [[6]]
## [1] "7] Sparsity = 0.530664365224743"
##
## [[7]]
## [1] "6] Any OTU sum to 1 or less? YES"
##
## [[8]]
## [1] "8] Number of singletons = 892"
##
## [[9]]
## [1] "9] Percent of OTUs that are singletons \n (i.e. exactly one read detected across all samples)0"
##
## [[10]]
## [1] "10] Number of sample variables are: 29"
##
## [[11]]
## [1] "Sample" "ID" "Depth"
## [4] "Depth_cm" "Site" "Sample_type"
## [7] "Date" "Location" "Temp...C..int"
## [10] "Salnidad..ppm.int" "Salinidad..ppt.int" "pH.int"
## [13] "pHsup" "Redox.medido.en.campo1" "Redox.mV..Eh.int"
## [16] "Temperatura...C.sup" "Salnidad..ppm.sup" "Salinidad..ppt.sup"
## [19] "Redox.medido.en.campo" "Redox.mV..Eh.sup" "S.2..mg.l..campo"
## [22] "ml.de.muestra1" "Factor.de.dilucion1" "S.2..mg.l."
## [25] "ID_M" "SO.....mg.l..campo" "ml.de.muestra2"
## [28] "Factor.de.dilucion2" "SO.....mg.l."
microbiome::summarize_phyloseq(ps_s6)## Compositional = NO2
## 1] Min. number of reads = 1549192] Max. number of reads = 2019623] Total number of reads = 10440904] Average number of reads = 1740155] Median number of reads = 172266.57] Sparsity = 0.656803344136026] Any OTU sum to 1 or less? YES8] Number of singletons = 29849] Percent of OTUs that are singletons
## (i.e. exactly one read detected across all samples)010] Number of sample variables are: 29SampleIDDepthDepth_cmSiteSample_typeDateLocationTemp...C..intSalnidad..ppm.intSalinidad..ppt.intpH.intpHsupRedox.medido.en.campo1Redox.mV..Eh.intTemperatura...C.supSalnidad..ppm.supSalinidad..ppt.supRedox.medido.en.campoRedox.mV..Eh.supS.2..mg.l..campoml.de.muestra1Factor.de.dilucion1S.2..mg.l.ID_MSO.....mg.l..campoml.de.muestra2Factor.de.dilucion2SO.....mg.l.2
## [[1]]
## [1] "1] Min. number of reads = 154919"
##
## [[2]]
## [1] "2] Max. number of reads = 201962"
##
## [[3]]
## [1] "3] Total number of reads = 1044090"
##
## [[4]]
## [1] "4] Average number of reads = 174015"
##
## [[5]]
## [1] "5] Median number of reads = 172266.5"
##
## [[6]]
## [1] "7] Sparsity = 0.65680334413602"
##
## [[7]]
## [1] "6] Any OTU sum to 1 or less? YES"
##
## [[8]]
## [1] "8] Number of singletons = 2984"
##
## [[9]]
## [1] "9] Percent of OTUs that are singletons \n (i.e. exactly one read detected across all samples)0"
##
## [[10]]
## [1] "10] Number of sample variables are: 29"
##
## [[11]]
## [1] "Sample" "ID" "Depth"
## [4] "Depth_cm" "Site" "Sample_type"
## [7] "Date" "Location" "Temp...C..int"
## [10] "Salnidad..ppm.int" "Salinidad..ppt.int" "pH.int"
## [13] "pHsup" "Redox.medido.en.campo1" "Redox.mV..Eh.int"
## [16] "Temperatura...C.sup" "Salnidad..ppm.sup" "Salinidad..ppt.sup"
## [19] "Redox.medido.en.campo" "Redox.mV..Eh.sup" "S.2..mg.l..campo"
## [22] "ml.de.muestra1" "Factor.de.dilucion1" "S.2..mg.l."
## [25] "ID_M" "SO.....mg.l..campo" "ml.de.muestra2"
## [28] "Factor.de.dilucion2" "SO.....mg.l."
microbiome::summarize_phyloseq(ps_s7)## Compositional = NO2
## 1] Min. number of reads = 664132] Max. number of reads = 2157403] Total number of reads = 12879034] Average number of reads = 160987.8755] Median number of reads = 174806.57] Sparsity = 0.597594053825566] Any OTU sum to 1 or less? YES8] Number of singletons = 22379] Percent of OTUs that are singletons
## (i.e. exactly one read detected across all samples)010] Number of sample variables are: 29SampleIDDepthDepth_cmSiteSample_typeDateLocationTemp...C..intSalnidad..ppm.intSalinidad..ppt.intpH.intpHsupRedox.medido.en.campo1Redox.mV..Eh.intTemperatura...C.supSalnidad..ppm.supSalinidad..ppt.supRedox.medido.en.campoRedox.mV..Eh.supS.2..mg.l..campoml.de.muestra1Factor.de.dilucion1S.2..mg.l.ID_MSO.....mg.l..campoml.de.muestra2Factor.de.dilucion2SO.....mg.l.2
## [[1]]
## [1] "1] Min. number of reads = 66413"
##
## [[2]]
## [1] "2] Max. number of reads = 215740"
##
## [[3]]
## [1] "3] Total number of reads = 1287903"
##
## [[4]]
## [1] "4] Average number of reads = 160987.875"
##
## [[5]]
## [1] "5] Median number of reads = 174806.5"
##
## [[6]]
## [1] "7] Sparsity = 0.59759405382556"
##
## [[7]]
## [1] "6] Any OTU sum to 1 or less? YES"
##
## [[8]]
## [1] "8] Number of singletons = 2237"
##
## [[9]]
## [1] "9] Percent of OTUs that are singletons \n (i.e. exactly one read detected across all samples)0"
##
## [[10]]
## [1] "10] Number of sample variables are: 29"
##
## [[11]]
## [1] "Sample" "ID" "Depth"
## [4] "Depth_cm" "Site" "Sample_type"
## [7] "Date" "Location" "Temp...C..int"
## [10] "Salnidad..ppm.int" "Salinidad..ppt.int" "pH.int"
## [13] "pHsup" "Redox.medido.en.campo1" "Redox.mV..Eh.int"
## [16] "Temperatura...C.sup" "Salnidad..ppm.sup" "Salinidad..ppt.sup"
## [19] "Redox.medido.en.campo" "Redox.mV..Eh.sup" "S.2..mg.l..campo"
## [22] "ml.de.muestra1" "Factor.de.dilucion1" "S.2..mg.l."
## [25] "ID_M" "SO.....mg.l..campo" "ml.de.muestra2"
## [28] "Factor.de.dilucion2" "SO.....mg.l."
Remove singletons
#Delete singletons
ps_s1 <- filter_taxa(ps_s1, function(x) sum(x) > 1, TRUE)
ps_s2 <- filter_taxa(ps_s2, function(x) sum(x) > 1, TRUE)
ps_s3 <- filter_taxa(ps_s3, function(x) sum(x) > 1, TRUE)
ps_s4 <- filter_taxa(ps_s4, function(x) sum(x) > 1, TRUE)
ps_s5 <- filter_taxa(ps_s5, function(x) sum(x) > 1, TRUE)
ps_s6 <- filter_taxa(ps_s6, function(x) sum(x) > 1, TRUE)
ps_s7 <- filter_taxa(ps_s7, function(x) sum(x) > 1, TRUE)
#check
ps_s1## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 6310 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 6310 taxa by 7 taxonomic ranks ]
ps_s2## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 5831 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 5831 taxa by 7 taxonomic ranks ]
ps_s3## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 4781 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 4781 taxa by 7 taxonomic ranks ]
ps_s4## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 5408 taxa and 4 samples ]
## sample_data() Sample Data: [ 4 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 5408 taxa by 7 taxonomic ranks ]
ps_s5## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 6205 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 6205 taxa by 7 taxonomic ranks ]
ps_s6## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 4113 taxa and 6 samples ]
## sample_data() Sample Data: [ 6 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 4113 taxa by 7 taxonomic ranks ]
ps_s7## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 4860 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 4860 taxa by 7 taxonomic ranks ]
library(microbiome)
library(microbiomeutilities)
#function
convert_to_compositional_and_index_reformat <- function(ps) {
ps_rel <- microbiome::transform(ps, "compositional")
ps_rel_f <- format_to_besthit(ps_rel)
return(ps_rel_f)
}
# phyloseq objects
phyloseq_objects <- list(ps_s1 = ps_s1, ps_s2 = ps_s2, ps_s3 = ps_s3,
ps_s4 = ps_s4, ps_s5 = ps_s5, ps_s6 = ps_s6,
ps_s7 = ps_s7, ps=ps)
# Apply
phyloseq_objects_rel <- lapply(phyloseq_objects, convert_to_compositional_and_index_reformat)
# Subset
ps_s1_rel_f <- phyloseq_objects_rel$ps_s1
ps_s2_rel_f <- phyloseq_objects_rel$ps_s2
ps_s3_rel_f <- phyloseq_objects_rel$ps_s3
ps_s4_rel_f <- phyloseq_objects_rel$ps_s4
ps_s5_rel_f <- phyloseq_objects_rel$ps_s5
ps_s6_rel_f <- phyloseq_objects_rel$ps_s6
ps_s7_rel_f <- phyloseq_objects_rel$ps_s7
ps_full_rel_f <- phyloseq_objects_rel$ps# Get full core with microbiome
full_ps_core_relative <- core(ps_full_rel_f, 0, 0.5)
full_ps_core_counts <- core(ps, 0, 0.5)
# Check point
ps## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 7097 taxa and 50 samples ]
## sample_data() Sample Data: [ 50 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 7097 taxa by 7 taxonomic ranks ]
ps_full_rel_f## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 7097 taxa and 50 samples ]
## sample_data() Sample Data: [ 50 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 7097 taxa by 8 taxonomic ranks ]
full_ps_core_relative## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 2048 taxa and 50 samples ]
## sample_data() Sample Data: [ 50 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 2048 taxa by 8 taxonomic ranks ]
full_ps_core_counts## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 2048 taxa and 50 samples ]
## sample_data() Sample Data: [ 50 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 2048 taxa by 7 taxonomic ranks ]
microbiome::summarize_phyloseq(full_ps_core_counts)## Compositional = NO2
## 1] Min. number of reads = 456802] Max. number of reads = 4514903] Total number of reads = 60921854] Average number of reads = 121843.75] Median number of reads = 1146517] Sparsity = 0.3219042968756] Any OTU sum to 1 or less? NO8] Number of singletons = 09] Percent of OTUs that are singletons
## (i.e. exactly one read detected across all samples)010] Number of sample variables are: 29SampleIDDepthDepth_cmSiteSample_typeDateLocationTemp...C..intSalnidad..ppm.intSalinidad..ppt.intpH.intpHsupRedox.medido.en.campo1Redox.mV..Eh.intTemperatura...C.supSalnidad..ppm.supSalinidad..ppt.supRedox.medido.en.campoRedox.mV..Eh.supS.2..mg.l..campoml.de.muestra1Factor.de.dilucion1S.2..mg.l.ID_MSO.....mg.l..campoml.de.muestra2Factor.de.dilucion2SO.....mg.l.2
## [[1]]
## [1] "1] Min. number of reads = 45680"
##
## [[2]]
## [1] "2] Max. number of reads = 451490"
##
## [[3]]
## [1] "3] Total number of reads = 6092185"
##
## [[4]]
## [1] "4] Average number of reads = 121843.7"
##
## [[5]]
## [1] "5] Median number of reads = 114651"
##
## [[6]]
## [1] "7] Sparsity = 0.321904296875"
##
## [[7]]
## [1] "6] Any OTU sum to 1 or less? NO"
##
## [[8]]
## [1] "8] Number of singletons = 0"
##
## [[9]]
## [1] "9] Percent of OTUs that are singletons \n (i.e. exactly one read detected across all samples)0"
##
## [[10]]
## [1] "10] Number of sample variables are: 29"
##
## [[11]]
## [1] "Sample" "ID" "Depth"
## [4] "Depth_cm" "Site" "Sample_type"
## [7] "Date" "Location" "Temp...C..int"
## [10] "Salnidad..ppm.int" "Salinidad..ppt.int" "pH.int"
## [13] "pHsup" "Redox.medido.en.campo1" "Redox.mV..Eh.int"
## [16] "Temperatura...C.sup" "Salnidad..ppm.sup" "Salinidad..ppt.sup"
## [19] "Redox.medido.en.campo" "Redox.mV..Eh.sup" "S.2..mg.l..campo"
## [22] "ml.de.muestra1" "Factor.de.dilucion1" "S.2..mg.l."
## [25] "ID_M" "SO.....mg.l..campo" "ml.de.muestra2"
## [28] "Factor.de.dilucion2" "SO.....mg.l."
Create full, core and Site ampvis objects
library(tibble)
library(ampvis2)
# Ampvis object function
create_ampvis_object <- function(pseq) {
# create and extract otu table
otu_table_ampvis <- data.frame(OTU = rownames(phyloseq::otu_table(pseq)@.Data),
phyloseq::otu_table(pseq)@.Data,
phyloseq::tax_table(pseq)@.Data,
check.names = FALSE)
# Metadata
meta_data_ampvis <- data.frame(phyloseq::sample_data(pseq),
check.names = FALSE)
# change index by SampleID
meta_data_ampvis <- meta_data_ampvis %>% rownames_to_column(var = "SampleID")
# ampvis object
av2 <- amp_load(otu_table_ampvis, meta_data_ampvis)
return(av2)
}# phyloseq objects
phyloseq_objects <-list(ps_s1 = ps_s1, ps_s2 = ps_s2, ps_s3 = ps_s3,
ps_s4 = ps_s4, ps_s5 = ps_s5, ps_s6 = ps_s6,
ps_s7 = ps_s7, ps=ps, full_ps_core_counts = full_ps_core_counts)
# Apply function
ampvis_objects <- lapply(phyloseq_objects, create_ampvis_object)
ampvis_objects## $ps_s1
## ampvis2 object with 3 elements.
## Summary of OTU table:
## Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads
## 8 6310 1429162 148635 218877 171227
## Avg#Reads
## 178645.25
##
## Assigned taxonomy:
## Kingdom Phylum Class Order Family Genus
## 6310(100%) 6291(99.7%) 6261(99.22%) 6055(95.96%) 6044(95.78%) 5946(94.23%)
## Species
## 4349(68.92%)
##
## Metadata variables: 30
## SampleID, Sample, ID, Depth, Depth_cm, Site, Sample_type, Date, Location, Temp...C..int, Salnidad..ppm.int, Salinidad..ppt.int, pH.int, pHsup, Redox.medido.en.campo1, Redox.mV..Eh.int, Temperatura...C.sup, Salnidad..ppm.sup, Salinidad..ppt.sup, Redox.medido.en.campo, Redox.mV..Eh.sup, S.2..mg.l..campo, ml.de.muestra1, Factor.de.dilucion1, S.2..mg.l., ID_M, SO.....mg.l..campo, ml.de.muestra2, Factor.de.dilucion2, SO.....mg.l.
##
## $ps_s2
## ampvis2 object with 3 elements.
## Summary of OTU table:
## Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads
## 8 5831 1379847 132348 220547 173469.5
## Avg#Reads
## 172480.88
##
## Assigned taxonomy:
## Kingdom Phylum Class Order Family Genus
## 5831(100%) 5812(99.67%) 5782(99.16%) 5581(95.71%) 5574(95.59%) 5481(94%)
## Species
## 4018(68.91%)
##
## Metadata variables: 30
## SampleID, Sample, ID, Depth, Depth_cm, Site, Sample_type, Date, Location, Temp...C..int, Salnidad..ppm.int, Salinidad..ppt.int, pH.int, pHsup, Redox.medido.en.campo1, Redox.mV..Eh.int, Temperatura...C.sup, Salnidad..ppm.sup, Salinidad..ppt.sup, Redox.medido.en.campo, Redox.mV..Eh.sup, S.2..mg.l..campo, ml.de.muestra1, Factor.de.dilucion1, S.2..mg.l., ID_M, SO.....mg.l..campo, ml.de.muestra2, Factor.de.dilucion2, SO.....mg.l.
##
## $ps_s3
## ampvis2 object with 3 elements.
## Summary of OTU table:
## Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads
## 8 4781 1686895 101979 674846 141534
## Avg#Reads
## 210861.88
##
## Assigned taxonomy:
## Kingdom Phylum Class Order Family Genus
## 4781(100%) 4764(99.64%) 4735(99.04%) 4557(95.31%) 4552(95.21%) 4477(93.64%)
## Species
## 3322(69.48%)
##
## Metadata variables: 30
## SampleID, Sample, ID, Depth, Depth_cm, Site, Sample_type, Date, Location, Temp...C..int, Salnidad..ppm.int, Salinidad..ppt.int, pH.int, pHsup, Redox.medido.en.campo1, Redox.mV..Eh.int, Temperatura...C.sup, Salnidad..ppm.sup, Salinidad..ppt.sup, Redox.medido.en.campo, Redox.mV..Eh.sup, S.2..mg.l..campo, ml.de.muestra1, Factor.de.dilucion1, S.2..mg.l., ID_M, SO.....mg.l..campo, ml.de.muestra2, Factor.de.dilucion2, SO.....mg.l.
##
## $ps_s4
## ampvis2 object with 3 elements.
## Summary of OTU table:
## Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads
## 4 5408 715609 127631 206504 190737
## Avg#Reads
## 178902.25
##
## Assigned taxonomy:
## Kingdom Phylum Class Order Family Genus
## 5408(100%) 5393(99.72%) 5366(99.22%) 5185(95.88%) 5176(95.71%) 5091(94.14%)
## Species
## 3694(68.31%)
##
## Metadata variables: 30
## SampleID, Sample, ID, Depth, Depth_cm, Site, Sample_type, Date, Location, Temp...C..int, Salnidad..ppm.int, Salinidad..ppt.int, pH.int, pHsup, Redox.medido.en.campo1, Redox.mV..Eh.int, Temperatura...C.sup, Salnidad..ppm.sup, Salinidad..ppt.sup, Redox.medido.en.campo, Redox.mV..Eh.sup, S.2..mg.l..campo, ml.de.muestra1, Factor.de.dilucion1, S.2..mg.l., ID_M, SO.....mg.l..campo, ml.de.muestra2, Factor.de.dilucion2, SO.....mg.l.
##
## $ps_s5
## ampvis2 object with 3 elements.
## Summary of OTU table:
## Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads
## 8 6205 1462992 97081 298225 167958
## Avg#Reads
## 182874
##
## Assigned taxonomy:
## Kingdom Phylum Class Order Family Genus
## 6205(100%) 6188(99.73%) 6159(99.26%) 5952(95.92%) 5943(95.78%) 5838(94.09%)
## Species
## 4243(68.38%)
##
## Metadata variables: 30
## SampleID, Sample, ID, Depth, Depth_cm, Site, Sample_type, Date, Location, Temp...C..int, Salnidad..ppm.int, Salinidad..ppt.int, pH.int, pHsup, Redox.medido.en.campo1, Redox.mV..Eh.int, Temperatura...C.sup, Salnidad..ppm.sup, Salinidad..ppt.sup, Redox.medido.en.campo, Redox.mV..Eh.sup, S.2..mg.l..campo, ml.de.muestra1, Factor.de.dilucion1, S.2..mg.l., ID_M, SO.....mg.l..campo, ml.de.muestra2, Factor.de.dilucion2, SO.....mg.l.
##
## $ps_s6
## ampvis2 object with 3 elements.
## Summary of OTU table:
## Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads
## 6 4113 1044090 154919 201962 172266.5
## Avg#Reads
## 174015
##
## Assigned taxonomy:
## Kingdom Phylum Class Order Family Genus
## 4113(100%) 4101(99.71%) 4084(99.29%) 3938(95.75%) 3927(95.48%) 3830(93.12%)
## Species
## 2702(65.69%)
##
## Metadata variables: 30
## SampleID, Sample, ID, Depth, Depth_cm, Site, Sample_type, Date, Location, Temp...C..int, Salnidad..ppm.int, Salinidad..ppt.int, pH.int, pHsup, Redox.medido.en.campo1, Redox.mV..Eh.int, Temperatura...C.sup, Salnidad..ppm.sup, Salinidad..ppt.sup, Redox.medido.en.campo, Redox.mV..Eh.sup, S.2..mg.l..campo, ml.de.muestra1, Factor.de.dilucion1, S.2..mg.l., ID_M, SO.....mg.l..campo, ml.de.muestra2, Factor.de.dilucion2, SO.....mg.l.
##
## $ps_s7
## ampvis2 object with 3 elements.
## Summary of OTU table:
## Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads
## 8 4860 1287903 66413 215740 174806.5
## Avg#Reads
## 160987.88
##
## Assigned taxonomy:
## Kingdom Phylum Class Order Family Genus
## 4860(100%) 4845(99.69%) 4826(99.3%) 4663(95.95%) 4651(95.7%) 4552(93.66%)
## Species
## 3230(66.46%)
##
## Metadata variables: 30
## SampleID, Sample, ID, Depth, Depth_cm, Site, Sample_type, Date, Location, Temp...C..int, Salnidad..ppm.int, Salinidad..ppt.int, pH.int, pHsup, Redox.medido.en.campo1, Redox.mV..Eh.int, Temperatura...C.sup, Salnidad..ppm.sup, Salinidad..ppt.sup, Redox.medido.en.campo, Redox.mV..Eh.sup, S.2..mg.l..campo, ml.de.muestra1, Factor.de.dilucion1, S.2..mg.l., ID_M, SO.....mg.l..campo, ml.de.muestra2, Factor.de.dilucion2, SO.....mg.l.
##
## $ps
## ampvis2 object with 3 elements.
## Summary of OTU table:
## Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads
## 50 7097 9006498 66413 674846 171486.5
## Avg#Reads
## 180129.96
##
## Assigned taxonomy:
## Kingdom Phylum Class Order Family Genus
## 7097(100%) 7077(99.72%) 7043(99.24%) 6811(95.97%) 6795(95.74%) 6672(94.01%)
## Species
## 4818(67.89%)
##
## Metadata variables: 30
## SampleID, Sample, ID, Depth, Depth_cm, Site, Sample_type, Date, Location, Temp...C..int, Salnidad..ppm.int, Salinidad..ppt.int, pH.int, pHsup, Redox.medido.en.campo1, Redox.mV..Eh.int, Temperatura...C.sup, Salnidad..ppm.sup, Salinidad..ppt.sup, Redox.medido.en.campo, Redox.mV..Eh.sup, S.2..mg.l..campo, ml.de.muestra1, Factor.de.dilucion1, S.2..mg.l., ID_M, SO.....mg.l..campo, ml.de.muestra2, Factor.de.dilucion2, SO.....mg.l.
##
## $full_ps_core_counts
## ampvis2 object with 3 elements.
## Summary of OTU table:
## Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads
## 50 2048 6092185 45680 451490 114651
## Avg#Reads
## 121843.7
##
## Assigned taxonomy:
## Kingdom Phylum Class Order Family Genus
## 2048(100%) 2043(99.76%) 2031(99.17%) 1960(95.7%) 1957(95.56%) 1913(93.41%)
## Species
## 1336(65.23%)
##
## Metadata variables: 30
## SampleID, Sample, ID, Depth, Depth_cm, Site, Sample_type, Date, Location, Temp...C..int, Salnidad..ppm.int, Salinidad..ppt.int, pH.int, pHsup, Redox.medido.en.campo1, Redox.mV..Eh.int, Temperatura...C.sup, Salnidad..ppm.sup, Salinidad..ppt.sup, Redox.medido.en.campo, Redox.mV..Eh.sup, S.2..mg.l..campo, ml.de.muestra1, Factor.de.dilucion1, S.2..mg.l., ID_M, SO.....mg.l..campo, ml.de.muestra2, Factor.de.dilucion2, SO.....mg.l.
#subset amp obj
list(ps_s1 = ps_s1, ps_s2 = ps_s2, ps_s3 = ps_s3,
ps_s4 = ps_s4, ps_s5 = ps_s5, ps_s6 = ps_s6,
ps_s7 = ps_s7, ps=ps)## $ps_s1
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 6310 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 6310 taxa by 7 taxonomic ranks ]
##
## $ps_s2
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 5831 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 5831 taxa by 7 taxonomic ranks ]
##
## $ps_s3
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 4781 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 4781 taxa by 7 taxonomic ranks ]
##
## $ps_s4
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 5408 taxa and 4 samples ]
## sample_data() Sample Data: [ 4 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 5408 taxa by 7 taxonomic ranks ]
##
## $ps_s5
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 6205 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 6205 taxa by 7 taxonomic ranks ]
##
## $ps_s6
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 4113 taxa and 6 samples ]
## sample_data() Sample Data: [ 6 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 4113 taxa by 7 taxonomic ranks ]
##
## $ps_s7
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 4860 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 4860 taxa by 7 taxonomic ranks ]
##
## $ps
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 7097 taxa and 50 samples ]
## sample_data() Sample Data: [ 50 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 7097 taxa by 7 taxonomic ranks ]
av2_full <- ampvis_objects$ps
av2_full_core <- ampvis_objects$full_ps_core_counts
av2_s1 <- ampvis_objects$ps_s1
av2_s2 <-ampvis_objects$ps_s2
av2_s3 <- ampvis_objects$ps_s3
av2_s4 <- ampvis_objects$ps_s4
av2_s5 <- ampvis_objects$ps_s5
av2_s6 <- ampvis_objects$ps_s6
av2_s7 <- ampvis_objects$ps_s7All
#Factor Depth_cm
av2_full$metadata$Depth <- factor(av2_full$metadata$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
# Full all samples
ampv_heatmap_abundances_phylum_full <- amp_heatmap(av2_full,
group_by = "Depth",
facet_by = "Site",
plot_values = TRUE,
tax_show = 25,
tax_aggregate = "Phylum",
tax_add = "Kingdom",
plot_colorscale = "log10",
plot_legendbreaks = c(1, 5, 5),
color_vector = c("white","deepskyblue3","magenta3"))
ampv_heatmap_abundances_phylum_fullCore
#Factor Depth_cm
av2_full_core$metadata$Depth <- factor(av2_full_core$metadata$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
# Full all samples
ampv_heatmap_abundances_phylum_core <- amp_heatmap(av2_full_core,
group_by = "Depth",
facet_by = "Site",
plot_values = TRUE,
tax_show = 25,
tax_aggregate = "Phylum",
tax_add = "Kingdom",
plot_colorscale = "log10",
plot_legendbreaks = c(1, 5, 5),
color_vector = c("white","deepskyblue3","magenta3"))
ampv_heatmap_abundances_phylum_coreCombine plots
library(cowplot)
#Combine venn plot
title_abund_plot <- ggdraw() + draw_label("Relative Abundance", fontface = 'bold', x = 0.5, hjust = 0.5, size = 15)
ampv_all_g_plots <- plot_grid(title_abund_plot, ampv_heatmap_abundances_phylum_full, ampv_heatmap_abundances_phylum_core, labels = c("", "A", "B"), ncol = 1, rel_heights = c(0.15, 1, 1))
ampv_all_g_plotsAll
library(ampvis2)
library(ggplot2)
#Factor Depth_cm
av2_full$metadata$Depth <- factor(av2_full$metadata$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
# rank abundance
rank_ab_full <- amp_rankabundance(av2_full, log10_x = TRUE, group_by = "Depth") +
scale_color_lancet()
rank_ab_full#ggsave("results/plots/02.diversity/02.rank_abundance.png", rank_ab , width = 5, height = 4)Site
# Site1
#Factor Depth_cm
av2_s1$metadata$Depth <- factor(av2_s1$metadata$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
# rank abundance
rank_ab_s1 <- amp_rankabundance(av2_s1, log10_x = TRUE, group_by = "Depth") +
scale_color_lancet()
# Site2
#Factor Depth_cm
av2_s2$metadata$Depth <- factor(av2_s2$metadata$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
# rank abundance
rank_ab_s2 <- amp_rankabundance(av2_s2, log10_x = TRUE, group_by = "Depth") +
scale_color_lancet()
# Site3
#Factor Depth_cm
av2_s3$metadata$Depth <- factor(av2_s3$metadata$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
# rank abundance
rank_ab_s3 <- amp_rankabundance(av2_s3, log10_x = TRUE, group_by = "Depth") +
scale_color_lancet()
# Site4
#Factor Depth_cm
av2_s4$metadata$Depth <- factor(av2_s4$metadata$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
# rank abundance
rank_ab_s4 <- amp_rankabundance(av2_s4, log10_x = TRUE, group_by = "Depth") +
scale_color_lancet()
# Site5
#Factor Depth_cm
av2_s5$metadata$Depth <- factor(av2_s5$metadata$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
# rank abundance
rank_ab_s5 <- amp_rankabundance(av2_s5, log10_x = TRUE, group_by = "Depth") +
scale_color_lancet()
# Site6
#Factor Depth_cm
av2_s6$metadata$Depth <- factor(av2_s6$metadata$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
# rank abundance
rank_ab_s6 <- amp_rankabundance(av2_s6, log10_x = TRUE, group_by = "Depth") +
scale_color_lancet()
# Site7
#Factor Depth_cm
av2_s7$metadata$Depth <- factor(av2_s7$metadata$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
# rank abundance
rank_ab_s7 <- amp_rankabundance(av2_s7, log10_x = TRUE, group_by = "Depth") +
scale_color_lancet()Core
library(ampvis2)
library(ggplot2)
#Factor Depth_cm
av2_full_core$metadata$Depth <- factor(av2_full_core$metadata$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
# rank abundance
rank_ab_full_core <- amp_rankabundance(av2_full_core, log10_x = TRUE, group_by = "Depth") +
scale_color_lancet()
rank_ab_full_coreCombine
library(cowplot)
#Combine euler plot
title_rank_plot <- ggdraw() + draw_label("Relative diversity and species evenness", fontface = 'bold', x = 0.5, hjust = 0.5, size = 15)
rank_ab_full_core_s1 <- plot_grid(rank_ab_full, rank_ab_full_core, rank_ab_s1, labels = c("F", "C", "S1"), ncol = 3)
rank_ab_s2tos4 <- plot_grid(rank_ab_s2,rank_ab_s3, rank_ab_s4, labels = c("S2","S3", "S4"), ncol = 3)## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
rank_ab_s5tos7 <- plot_grid(rank_ab_s5, rank_ab_s6, rank_ab_s7, labels = c("S5", "S6", "S7"), ncol = 3)
rankabund_all_plots <- plot_grid(title_rank_plot, rank_ab_full_core_s1, rank_ab_s2tos4, rank_ab_s5tos7, ncol = 1, rel_heights = c(0.14, 1, 1, 1))
rankabund_all_plotslibrary(hilldiv)## Registered S3 methods overwritten by 'FSA':
## method from
## confint.boot car
## hist.boot car
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ readr 2.1.5
## ✔ lubridate 1.9.3 ✔ stringr 1.5.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ microbiome::alpha() masks ggplot2::alpha()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ lubridate::stamp() masks cowplot::stamp()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(kableExtra)##
## Attaching package: 'kableExtra'
##
## The following object is masked from 'package:dplyr':
##
## group_rows
# Get hill numbers
q0 <- hill_div(otu_table, qvalue = 0)
q1 <- hill_div(otu_table, qvalue = 1)
q2 <- hill_div(otu_table, qvalue = 2)
# Merge metadata with Hill numbers
q012_all <- cbind(q0, q1, q2) %>% as.data.frame() %>% rownames_to_column(var = "SampleID")
metadata_with_hill <- q012_all %>%
inner_join(metadata, by = c("SampleID"="SampleID"))
#save table
write.table(metadata_with_hill, "../results/05.diversity/Metadata_with_hill.tsv", quote = FALSE, sep = "\t", row.names = TRUE, col.names=TRUE)
# show
metadata_with_hill %>% head() %>% kable()| SampleID | q0 | q1 | q2 | Sample | ID | Depth | Depth_cm | Site | Sample_type | Date | Location | Temp (°C) int | Salnidad (ppm)int | Salinidad (ppt)int | pH int | pHsup | Redox medido en campo1 | Redox mV (Eh)int | Temperatura (°C)sup | Salnidad (ppm)sup | Salinidad (ppt)sup | Redox medido en campo | Redox mV (Eh)sup | S-2 (mg/l) campo | ml de muestra1 | Factor de dilucion1 | S-2 (mg/l) | ID_M | SO₄²- (mg/l) campo | ml de muestra2 | Factor de dilucion2 | SO₄ ² (mg/l) |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| S1A015 | 2267 | 581.0228 | 172.50572 | S1A015 | S1A | 1 | 0-15 | Site_1 | Sediment | 28/09/2023 | Dique Miguelito | 31.4 | 817.5 | 0.82 | 6.55 | 7.34 | -114 | 121 | 33.0 | 734.7 | 0.73 | -108 | 127 | 7.73 | 25 | 1 | 7.73 | S1a M | 53.7 | 2 | 10 | 537 |
| S1A1530 | 3082 | 924.7261 | 341.75340 | S1A1530 | S1A | 2 | 16-30 | Site_1 | Sediment | 28/09/2023 | Dique Miguelito | 31.4 | 817.5 | 0.82 | 6.55 | 7.34 | -114 | 121 | 33.0 | 734.7 | 0.73 | -108 | 127 | 7.73 | 25 | 1 | 7.73 | S1a M | 53.7 | 2 | 10 | 537 |
| S1A3050 | 3234 | 926.0288 | 303.32359 | S1A3050 | S1A | 3 | 31-45 | Site_1 | Sediment | 28/09/2023 | Dique Miguelito | 31.4 | 817.5 | 0.82 | 6.55 | 7.34 | -114 | 121 | 33.0 | 734.7 | 0.73 | -108 | 127 | 7.73 | 25 | 1 | 7.73 | S1a M | 53.7 | 2 | 10 | 537 |
| S1A5075 | 3506 | 879.5831 | 228.44914 | S1A5075 | S1A | 4 | 50-75 | Site_1 | Sediment | 28/09/2023 | Dique Miguelito | 31.4 | 817.5 | 0.82 | 6.55 | 7.34 | -114 | 121 | 33.0 | 734.7 | 0.73 | -108 | 127 | 7.73 | 25 | 1 | 7.73 | S1a M | 53.7 | 2 | 10 | 537 |
| S1B07 | 3011 | 364.4166 | 62.98918 | S1B07 | S1B | 1 | 0-15 | Site_1 | Sediment | 28/09/2023 | Dique Miguelito | 30.9 | 794.0 | 0.79 | 6.69 | 7.79 | -234 | 1 | 33.2 | 732.8 | 0.73 | -136 | 99 | 7.09 | 25 | 1 | 7.09 | S1b M | 55.0 | 2 | 10 | 550 |
| S1B2230 | 3973 | 1127.1162 | 320.30757 | S1B2230 | S1B | 2 | 16-30 | Site_1 | Sediment | 28/09/2023 | Dique Miguelito | 30.9 | 794.0 | 0.79 | 6.69 | 7.79 | -234 | 1 | 33.2 | 732.8 | 0.73 | -136 | 99 | 7.09 | 25 | 1 | 7.09 | S1b M | 55.0 | 2 | 10 | 550 |
Get means
# Reorder q values
meta_qs <- metadata_with_hill %>%
pivot_longer(cols = q0:q2, names_to = "q", values_to = "value") %>%
filter(q %in% c("q0", "q1", "q2")) %>%
mutate(
qs = case_when(
q == "q0" ~ "q=0",
q == "q1" ~ "q=1",
q == "q2" ~ "q=2"
))
#Get means of Hill numbers
means <- meta_qs %>% group_by(Depth_cm, Site, qs) %>%
summarise(means = mean(value, na.rm = TRUE),
sd = sd(value, na.rm = TRUE),
.groups = 'drop')
#save table
write.table(means, "../results/05.diversity/Hill_means_sd.tsv", quote = FALSE, sep = "\t", row.names = TRUE, col.names=TRUE)
print(means)## # A tibble: 81 × 5
## Depth_cm Site qs means sd
## <chr> <chr> <chr> <dbl> <dbl>
## 1 0-15 Site_1 q=0 2639 526.
## 2 0-15 Site_1 q=1 473. 153.
## 3 0-15 Site_1 q=2 118. 77.4
## 4 0-15 Site_2 q=0 3142. 771.
## 5 0-15 Site_2 q=1 941. 339.
## 6 0-15 Site_2 q=2 288. 116.
## 7 0-15 Site_3 q=0 2771 781.
## 8 0-15 Site_3 q=1 612. 429.
## 9 0-15 Site_3 q=2 150. 123.
## 10 0-15 Site_4 q=0 3009 NA
## # ℹ 71 more rows
library(ggpubr) ##
## Attaching package: 'ggpubr'
## The following object is masked from 'package:cowplot':
##
## get_legend
## The following object is masked from 'package:qiime2R':
##
## mean_sd
# factor to reorder plot
metadata_with_hill$Site <- factor(metadata_with_hill$Site)
metadata_with_hill$Depth <- factor(metadata_with_hill$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
#plot
hill_barplot_explore <- metadata_with_hill %>%
pivot_longer(cols = q0:q2, names_to = "q", values_to = "value") %>%
filter(q %in% c("q0", "q1", "q2")) %>%
mutate(
qs = case_when(
q == "q0" ~ "q=0",
q == "q1" ~ "q=1",
q == "q2" ~ "q=2"
)) %>%
ggbarplot(
x = "Depth",
y = "value",
add = "mean_se",
facet.by = c("qs", "Site"),
fill = "Depth"
) +
scale_fill_lancet() +
#geom_jitter(size = 0.8, position = position_jitter(width = 0.1)) +
facet_grid(rows = vars(qs), cols = vars(Site), scales = "free_y") +
theme(
strip.background = element_blank(), #element_rect(fill = "")
strip.text.x = element_text(face= "bold", size = 13, margin = margin(0.5, 0, 0.5, 0)),
strip.text.y = element_text(face= "bold", size = 14, angle = 0, margin = margin(0, 0.5, 0, 0.5)),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.line = element_line(colour = "black"),
#panel.border = element_blank(),
panel.spacing.x = unit(0.5, "lines"),
panel.background = element_blank(),
axis.title.y = element_blank(),
axis.text.y = element_text(colour = "black", size = 11),
legend.position = "bottom"
) +
labs(
fill = "Depth",
y = "",
x = "",
title = "Biodiversity across Sites and Depths"
) +
theme(legend.text = element_text(size = 12))## Warning in stats::qt(ci/2 + 0.5, data_sum$length - 1): NaNs produced
## Warning in stats::qt(ci/2 + 0.5, data_sum$length - 1): NaNs produced
## Warning in stats::qt(ci/2 + 0.5, data_sum$length - 1): NaNs produced
## Warning in stats::qt(ci/2 + 0.5, data_sum$length - 1): NaNs produced
## Warning in stats::qt(ci/2 + 0.5, data_sum$length - 1): NaNs produced
## Warning in stats::qt(ci/2 + 0.5, data_sum$length - 1): NaNs produced
## Warning in stats::qt(ci/2 + 0.5, data_sum$length - 1): NaNs produced
## Warning in stats::qt(ci/2 + 0.5, data_sum$length - 1): NaNs produced
## Warning in stats::qt(ci/2 + 0.5, data_sum$length - 1): NaNs produced
## Warning in stats::qt(ci/2 + 0.5, data_sum$length - 1): NaNs produced
## Warning in stats::qt(ci/2 + 0.5, data_sum$length - 1): NaNs produced
## Warning in stats::qt(ci/2 + 0.5, data_sum$length - 1): NaNs produced
hill_barplot_explore#ggsave("../results/plots/Hill_diversity_across_Sites_and_Depths_fltr.png", hill_barplot_explore_fltr, width = 11, height = 11)Alpha diversity depth correlation to order q=0
#Get effective reads
Effective_reads <- colSums(otu_table[-1, ])
metadata_with_hill$Effective_reads <- Effective_reads#install.packages("ggpubr")
library(ggpubr)
q0_vs_depth <- ggscatter(metadata_with_hill,
x = "Effective_reads",
y = "q0",
xlab= "Sequencing depth (number of reads)",
ylab = "q=0",
#ylab="Alpha diversity q=0 (effective number of total ASVs)",
add = "reg.line", # Add regression line
add.params = list(color = "#114e9d", fill = "lightgray"), # Customize reg. line
conf.int = TRUE, # Add confidence interval
cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
cor.coeff.args = list(method = "pearson", label.x = 3, label.sep = "\n")
)+
theme(legend.title = element_blank(), legend.position = "none")Alpha diversity depth correlation to order q=1
q1_vs_depth <- ggscatter(metadata_with_hill,
x = "Effective_reads",
y = "q1",
xlab= "Sequencing depth (number of reads)",
ylab = "q=1",
#ylab="Alpha diversity q=1 (effective number of total ASVs)",
add = "reg.line", # Add regression line
add.params = list(color = "#114e9d", fill = "lightgray"), # Customize reg. line
conf.int = TRUE, # Add confidence interval
cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
cor.coeff.args = list(method = "pearson", label.x = 3, label.sep = "\n")
)+
theme(legend.title = element_blank(), legend.position = "none")Alpha diversity depth correlation to order q=2
q2_vs_depth <- ggscatter(metadata_with_hill,
x = "Effective_reads",
y = "q2",
xlab= "Sequencing depth (number of reads)",
ylab = "q=2",
#ylab="Alpha diversity q=2 (effective number of total ASVs)",
add = "reg.line", # Add regression line
add.params = list(color = "#114e9d", fill = "lightgray"), # Customize reg. line
conf.int = TRUE, # Add confidence interval
cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
cor.coeff.args = list(method = "pearson", label.x = 3, label.sep = "\n")
)+
theme(legend.title = element_blank(), legend.position = "none")library(ggplot2)
library(cowplot)
#save plots
#Combine plot
title_corr_plot <- ggdraw() + draw_label("Alpha diversity depth correlation with fltr samples")
correlation_plot_q012 <- plot_grid(title_corr_plot, q0_vs_depth, q1_vs_depth, q2_vs_depth, labels = c(" ","A", "B", "C"), ncol = 1, rel_heights = c(0.15, 1, 1, 1))
correlation_plot_q012#save plot
#ggsave("../results/plots/alpha_depth_correlation_samples.png", correlation_plot_q012_, width = 8, height = 8)Shapiro test to Hill numbers
library(ggplot2)
library(cowplot)
#Shapiro test
shapiro_q0 <- shapiro.test(metadata_with_hill$q0)
shapiro_q0_pvalue <- round(shapiro_q0$p.value, 5)
shapiro_q1 <- shapiro.test(metadata_with_hill$q1)
shapiro_q1_pvalue <- round(shapiro_q1$p.value, 5)
shapiro_q2 <- shapiro.test(metadata_with_hill$q2)
shapiro_q2_pvalue <- round(shapiro_q2$p.value, 5)
#Histograms
histplot_q0 <- ggplot(metadata_with_hill, aes(x = q0, xlab="q=0")) +
geom_histogram(fill = "lightblue", bins = 15) +
ggtitle(paste("Shapiro, p-value:", shapiro_q0_pvalue)) +
theme_bw() + xlab("q=0") + ylab("Frequency")
histplot_q1 <- ggplot(metadata_with_hill, aes(x = q1)) +
geom_histogram(fill = "lightblue", bins = 15) +
ggtitle(paste("Shapiro, p-value:", shapiro_q1_pvalue)) +
theme_bw() + xlab("q=1") + ylab("Frequency")
histplot_q2 <- ggplot(metadata_with_hill, aes(x = q2)) +
geom_histogram(fill = "lightblue", bins = 15) +
ggtitle(paste("Shapiro, p-value:", shapiro_q2_pvalue)) +
theme_bw() + xlab("q=2") + ylab("Frequency")
#Combine plot
title_plot <- ggdraw() + draw_label("Histogram of Hill diversity") #, fontface = 'bold', x = 0.5, hjust = 0.5)
histplot_q012 <- plot_grid(title_plot, histplot_q0, histplot_q1, histplot_q2, labels = c(" ","A", "B", "C"), ncol = 1, rel_heights = c(0.15, 1, 1, 1))
histplot_q012#save plot
#ggsave("../results/plots/hill_normality_test_histogram_samples.png", histplot_q012)A pesar de que tienen una distribución normal, dado que los datos del proyecto tienen medidas repetidas y desbalanceo, lo adecuado es hacer comparaciones basadas en modelos. En este caso modelos lineales mixtos
Estos modelos son útiles cuando las observaciones tienen estructuras de dependencia, como datos longitudinales o datos agrupados/anidados. Los modelos lineales mixtos incluyen tanto efectos fijos como efectos aleatorios, donde los efectos fijos son comunes a todas las unidades/individuos, y los efectos aleatorios varÃan por grupo o unidad.
lme
library(nlme)##
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
##
## collapse
library(pgirmess)## Registered S3 method overwritten by 'spdep':
## method from
## plot.mst ape
# lme
# q0
q0_lme <- lme(q0 ~ log(Effective_reads)+Depth_cm+Site, random=~1 |ID, data = metadata_with_hill)
q0_lme_perm <- PermTest(q0_lme, B = 100)
# q1
q1_lme <- lme(q1 ~ log(Effective_reads)+Depth_cm+Location, random=~1 |ID, data = metadata_with_hill)
q1_lme_perm <- PermTest(q1_lme, B = 100)
# q2
q2_lme <- lme(q2 ~ log(Effective_reads)+Depth_cm+Location, random=~1 |ID, data = metadata_with_hill)
q2_lme_perm <- PermTest(q2_lme, B = 100)Effect
# extract p-values fun
extract_pvalues <- function(perm_test_result) {
perm_test_result$resultats$p.value
}
# extract
lme_pvalues_q0 <- extract_pvalues(q0_lme_perm)
lme_pvalues_q1 <- extract_pvalues(q1_lme_perm)
lme_pvalues_q2 <- extract_pvalues(q2_lme_perm)
# Dataframe p-values
library(dplyr)
lme_pvalues_table <- data.frame(
Hill = c("q0", "q1", "q2"),
`Intercept` = c(lme_pvalues_q0[1], lme_pvalues_q1[1], lme_pvalues_q2[1]),
`log(Depth)` = c(lme_pvalues_q0[2], lme_pvalues_q1[2], lme_pvalues_q2[2]),
`Depth` = c(lme_pvalues_q0[3], lme_pvalues_q1[3], lme_pvalues_q2[3]),
`Site` = c(lme_pvalues_q0[4], lme_pvalues_q1[4], lme_pvalues_q2[4])
)
# table
knitr::kable(t(lme_pvalues_table), caption = "P-values from Permuted LME")| Hill | q0 | q1 | q2 |
| Intercept | 0.17 | 0.63 | 0.61 |
| log.Depth. | 0.00 | 0.95 | 0.81 |
| Depth | 0.00 | 0.14 | 0.08 |
| Site | 0.20 | 0.81 | 0.99 |
# save
write.csv(t(lme_pvalues_table), "../results/05.diversity/P-values_from_Permuted_LME.csv", row.names = FALSE)q0_lme_perm##
## Monte-Carlo test
##
## Call:
## PermTest.lme(obj = q0_lme, B = 100)
##
## Based on 100 replicates
## Simulated p-value:
## p.value
## (Intercept) 0.17
## log(Effective_reads) 0.00
## Depth_cm 0.00
## Site 0.20
Check with easystats
q=0
library(easystats)## # Attaching packages: easystats 0.7.1 (red = needs update)
## ✔ bayestestR 0.13.2 ✔ correlation 0.8.4
## ✔ datawizard 0.10.0 ✔ effectsize 0.8.7
## ✔ insight 0.19.10 ✔ modelbased 0.8.7
## ✔ performance 0.11.0 ✔ parameters 0.21.6
## ✔ report 0.5.8 ✖ see 0.8.3
##
## Restart the R-Session and update packages with `easystats::easystats_update()`.
check_q0_lme_1 <- check_model(q0_lme)## Converting missing values (`NA`) into regular values currently not
## possible for variables of class `NULL`.
print(check_q0_lme_1)q=1
check_q1_lme_1 <- check_model(q1_lme)## Converting missing values (`NA`) into regular values currently not
## possible for variables of class `NULL`.
print(check_q1_lme_1)q=2
check_q2_lme_1 <- check_model(q2_lme)## Converting missing values (`NA`) into regular values currently not
## possible for variables of class `NULL`.
check_q2_lme_1glm quasipoisson
Prepare data
aldex_clr_data <- t(getMonteCarloSample(aldex_clr,1))
pca <- prcomp(aldex_clr_data)
meta=metadata_with_hill %>% dplyr::rename(SampleID="SampleID")PCA
library(ggplot2)
library(dplyr)
library(tidyr)
# Reorder
meta$Depth <- factor(meta$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
# plot
pca_plot <- ggplot(data = data.frame(pca$x) %>%
rownames_to_column(var = "SampleID") %>%
left_join(meta, by = "SampleID")) +
geom_segment(data = data.frame(pca$rotation) %>%
rownames_to_column(var = "Feature.ID") %>%
mutate(a = sqrt(PC1^2 + PC2^2)) %>%
top_n(8, a) %>%
mutate(PC1 = PC1*50, PC2 = PC2*50),
aes(x = 0, xend = PC1, y = 0, yend = PC2)) +
geom_point(aes(x = PC1, y = PC2, color = Depth, shape = Site), size = 3) +
geom_vline(xintercept = 0, linetype = 2) +
geom_hline(yintercept = 0, linetype = 2) +
scale_color_lancet() +
scale_shape_manual(values = c(15, 16, 17, 18, 23, 24, 25)) +
scale_fill_lancet() +
theme(axis.text = element_text(colour = "black", size = 12),
axis.title = element_text(colour = "black", size = 12),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right",
legend.box = "vertical") +
labs(y = "PC2", x = "PC1", color = "Depth", shape = "Site") +
theme_bw()
pca_plot# Bray NMDS bray
library(vegan)## Loading required package: permute
##
## Attaching package: 'permute'
## The following object is masked from 'package:devtools':
##
## check
## Loading required package: lattice
## This is vegan 2.6-4
##
## Attaching package: 'vegan'
## The following object is masked from 'package:microbiome':
##
## diversity
bray=vegdist(t(otu_table), method = "bray")
nmds_source_bray = vegan::metaMDS(bray, trymax = 20, k = 2) ## Run 0 stress 0.1712067
## Run 1 stress 0.1756919
## Run 2 stress 0.1712066
## ... New best solution
## ... Procrustes: rmse 0.0001310485 max resid 0.0008477333
## ... Similar to previous best
## Run 3 stress 0.1742446
## Run 4 stress 0.1850122
## Run 5 stress 0.1885669
## Run 6 stress 0.1712066
## ... Procrustes: rmse 0.00003115916 max resid 0.0001995856
## ... Similar to previous best
## Run 7 stress 0.1712066
## ... Procrustes: rmse 0.00001631054 max resid 0.0001018519
## ... Similar to previous best
## Run 8 stress 0.1946932
## Run 9 stress 0.1742446
## Run 10 stress 0.1756918
## Run 11 stress 0.1756918
## Run 12 stress 0.1712068
## ... Procrustes: rmse 0.0002026832 max resid 0.001330524
## ... Similar to previous best
## Run 13 stress 0.1874831
## Run 14 stress 0.1742447
## Run 15 stress 0.1850122
## Run 16 stress 0.1900205
## Run 17 stress 0.1757299
## Run 18 stress 0.1741713
## Run 19 stress 0.1777425
## Run 20 stress 0.1712067
## ... Procrustes: rmse 0.000111613 max resid 0.0007341006
## ... Similar to previous best
## *** Best solution repeated 5 times
var_stress_nmds_bray <- round(nmds_source_bray$stress, 5)
var_stress_nmds_bray## [1] 0.17121
scores_source_bray= nmds_source_bray %>% vegan::scores()
metadata$Depth <- factor(metadata$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
nmds_plot_bray<- ggplot() +
geom_point(data=data.frame(scores_source_bray) %>%
rownames_to_column(var = "SampleID")%>%
left_join(metadata, by = "SampleID"),
aes(x=NMDS1, y=NMDS2, color = Depth, shape = Site), size=5) +
theme_linedraw()+
scale_fill_lancet()+
scale_color_lancet()+
scale_shape_manual(values = c(15, 16, 17, 18, 23, 24, 25)) +
theme(axis.text = element_text(colour = "black", size = 12),
axis.title = element_text(colour = "black", size = 12),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right",
legend.box = "vertical")+ theme_bw() +
labs(title="NMDS, Bray-Curtis distances") +
ylab("NMDS2")+xlab("NMDS1")
nmds_plot_bray <- nmds_plot_bray +
annotate("text", x = Inf, y = -Inf, label = paste("Stress:", var_stress_nmds_bray),
hjust = 3.1, vjust = -1.9, size = 3)
nmds_plot_bray# NMDS aichitson
library(vegan)
pseudocount <- 1
otu_table_pseudo <- otu_table + pseudocount
aitch=vegdist(t(otu_table_pseudo), method = "aitchison")
nmds_source_aitch = vegan::metaMDS(aitch, trymax = 20, k = 2) ## Run 0 stress 0.1034177
## Run 1 stress 0.1034177
## ... Procrustes: rmse 0.00000122641 max resid 0.000003746513
## ... Similar to previous best
## Run 2 stress 0.1542678
## Run 3 stress 0.1519305
## Run 4 stress 0.1034177
## ... Procrustes: rmse 0.000005404249 max resid 0.00002283729
## ... Similar to previous best
## Run 5 stress 0.1034177
## ... Procrustes: rmse 0.000001584979 max resid 0.000005991313
## ... Similar to previous best
## Run 6 stress 0.1519474
## Run 7 stress 0.1056559
## Run 8 stress 0.1056635
## Run 9 stress 0.1080944
## Run 10 stress 0.1034177
## ... Procrustes: rmse 0.000002979132 max resid 0.00001396478
## ... Similar to previous best
## Run 11 stress 0.1034177
## ... New best solution
## ... Procrustes: rmse 0.000001200118 max resid 0.000003370918
## ... Similar to previous best
## Run 12 stress 0.1034177
## ... Procrustes: rmse 0.000001639774 max resid 0.000004650393
## ... Similar to previous best
## Run 13 stress 0.1080944
## Run 14 stress 0.1034177
## ... Procrustes: rmse 0.000001682944 max resid 0.000006160126
## ... Similar to previous best
## Run 15 stress 0.1056635
## Run 16 stress 0.1080944
## Run 17 stress 0.1034177
## ... New best solution
## ... Procrustes: rmse 0.00000116747 max resid 0.000003465607
## ... Similar to previous best
## Run 18 stress 0.1034177
## ... Procrustes: rmse 0.000001206626 max resid 0.000002998858
## ... Similar to previous best
## Run 19 stress 0.1080944
## Run 20 stress 0.1034177
## ... Procrustes: rmse 0.000001326051 max resid 0.000004422972
## ... Similar to previous best
## *** Best solution repeated 3 times
var_stress_nmds_aitch <- round(nmds_source_aitch$stress, 5)
var_stress_nmds_aitch## [1] 0.10342
scores_source_aitch= nmds_source_aitch %>% vegan::scores()
metadata$Depth <- factor(metadata$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
nmds_plot_aitch<- ggplot() +
geom_point(data=data.frame(scores_source_aitch) %>%
rownames_to_column(var = "SampleID")%>%
left_join(metadata, by = "SampleID"),
aes(x=NMDS1, y=NMDS2, color = Depth, shape = Site), size=5) +
# theme_linedraw()+
scale_fill_lancet()+
scale_color_lancet()+
scale_shape_manual(values = c(15, 16, 17, 18, 23, 24, 25)) +
theme(axis.text = element_text(colour = "black", size = 12),
axis.title = element_text(colour = "black", size = 12),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right",
legend.box = "vertical")+ theme_bw() +
labs(title="NMDS, Aitchison distances") +
ylab("NMDS2")+xlab("NMDS1")
nmds_plot_aitch <- nmds_plot_aitch +
annotate("text", x = Inf, y = -Inf, label = paste("Stress:", var_stress_nmds_aitch),
hjust = 3.1, vjust = -1.9, size = 3)
nmds_plot_aitchlibrary(vegan)
bray_dist =vegdist(t(otu_table), method = "bray")
pcoas=wcmdscale(bray_dist, eig = T)
metadata$Depth <- factor(metadata$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
pcoa_plot<- ggplot() +
geom_point(data=data.frame(pcoas$points) %>% #individuals
rownames_to_column(var = "SampleID")%>%
left_join(metadata, by = "SampleID"),
aes(x=Dim1, y=Dim2,color = Depth, shape = Site ), size=4) +
geom_vline(xintercept = 0, linetype = 2) + #lines-cross
geom_hline(yintercept = 0, linetype = 2) +
theme_linedraw()+
scale_fill_lancet()+
scale_color_lancet()+
scale_shape_manual(values = c(15, 16, 17, 18, 23, 24, 25)) +
theme(axis.text = element_text(colour = "black", size = 12),
axis.title = element_text(colour = "black", size = 12),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right",
legend.box = "vertical")+ theme_bw()+
labs(title="PCoA, Bray-Curtis distances")+
ylab("PCoA2")+xlab("PCoA1")
pcoa_plot#Combine plot
nmds_bray_aitch <- plot_grid(nmds_plot_bray, nmds_plot_aitch, labels = c("A", "B"), ncol = 2)
nmds_bray_aitchmeta_just = data.frame(t(otu_table_pseudo), check.names = F) %>% rownames_to_column(var = "SampleID") %>% inner_join(metadata_with_hill)## Joining with `by = join_by(SampleID)`
permanovas = adonis2(t(otu_table_pseudo) ~ log(Effective_reads)+Depth_cm*Site,
data=meta_just, method = "aitchison",
permutations = 999, strata = meta_just$ID)
permanovasPrepare data
Explore
microbiome::summarize_phyloseq(full_ps_core_counts)## Compositional = NO2
## 1] Min. number of reads = 456802] Max. number of reads = 4514903] Total number of reads = 60921854] Average number of reads = 121843.75] Median number of reads = 1146517] Sparsity = 0.3219042968756] Any OTU sum to 1 or less? NO8] Number of singletons = 09] Percent of OTUs that are singletons
## (i.e. exactly one read detected across all samples)010] Number of sample variables are: 29SampleIDDepthDepth_cmSiteSample_typeDateLocationTemp...C..intSalnidad..ppm.intSalinidad..ppt.intpH.intpHsupRedox.medido.en.campo1Redox.mV..Eh.intTemperatura...C.supSalnidad..ppm.supSalinidad..ppt.supRedox.medido.en.campoRedox.mV..Eh.supS.2..mg.l..campoml.de.muestra1Factor.de.dilucion1S.2..mg.l.ID_MSO.....mg.l..campoml.de.muestra2Factor.de.dilucion2SO.....mg.l.2
## [[1]]
## [1] "1] Min. number of reads = 45680"
##
## [[2]]
## [1] "2] Max. number of reads = 451490"
##
## [[3]]
## [1] "3] Total number of reads = 6092185"
##
## [[4]]
## [1] "4] Average number of reads = 121843.7"
##
## [[5]]
## [1] "5] Median number of reads = 114651"
##
## [[6]]
## [1] "7] Sparsity = 0.321904296875"
##
## [[7]]
## [1] "6] Any OTU sum to 1 or less? NO"
##
## [[8]]
## [1] "8] Number of singletons = 0"
##
## [[9]]
## [1] "9] Percent of OTUs that are singletons \n (i.e. exactly one read detected across all samples)0"
##
## [[10]]
## [1] "10] Number of sample variables are: 29"
##
## [[11]]
## [1] "Sample" "ID" "Depth"
## [4] "Depth_cm" "Site" "Sample_type"
## [7] "Date" "Location" "Temp...C..int"
## [10] "Salnidad..ppm.int" "Salinidad..ppt.int" "pH.int"
## [13] "pHsup" "Redox.medido.en.campo1" "Redox.mV..Eh.int"
## [16] "Temperatura...C.sup" "Salnidad..ppm.sup" "Salinidad..ppt.sup"
## [19] "Redox.medido.en.campo" "Redox.mV..Eh.sup" "S.2..mg.l..campo"
## [22] "ml.de.muestra1" "Factor.de.dilucion1" "S.2..mg.l."
## [25] "ID_M" "SO.....mg.l..campo" "ml.de.muestra2"
## [28] "Factor.de.dilucion2" "SO.....mg.l."
#Get otu and metadata table of core phyloseq object
# extract otu table
otu_table_full_core <- otu_table(full_ps_core_counts@otu_table)
meta_full_core <- sample_data(full_ps_core_counts)# NMDS aichitson
library(vegan)
pseudocount <- 1
otu_table_full_core_pseudo <- otu_table_full_core + pseudocount
aitch_fc=vegdist(t(otu_table_full_core_pseudo), method = "aitchison")
nmds_source_aitch_fc = vegan::metaMDS(aitch_fc, trymax = 20, k = 2) ## Run 0 stress 0.09967979
## Run 1 stress 0.09967979
## ... Procrustes: rmse 0.000007858841 max resid 0.0000363496
## ... Similar to previous best
## Run 2 stress 0.1509729
## Run 3 stress 0.1023438
## Run 4 stress 0.09967979
## ... New best solution
## ... Procrustes: rmse 0.000004178664 max resid 0.0000199997
## ... Similar to previous best
## Run 5 stress 0.09967979
## ... Procrustes: rmse 0.000005413533 max resid 0.00002386895
## ... Similar to previous best
## Run 6 stress 0.1023717
## Run 7 stress 0.1023438
## Run 8 stress 0.1023438
## Run 9 stress 0.1722353
## Run 10 stress 0.09967979
## ... Procrustes: rmse 0.000001449633 max resid 0.000004770371
## ... Similar to previous best
## Run 11 stress 0.09967979
## ... Procrustes: rmse 0.000004391546 max resid 0.0000213067
## ... Similar to previous best
## Run 12 stress 0.09967979
## ... Procrustes: rmse 0.000002657153 max resid 0.00001148891
## ... Similar to previous best
## Run 13 stress 0.1023438
## Run 14 stress 0.09967979
## ... Procrustes: rmse 0.000001907427 max resid 0.000007322214
## ... Similar to previous best
## Run 15 stress 0.09967979
## ... Procrustes: rmse 0.000006794376 max resid 0.00003191459
## ... Similar to previous best
## Run 16 stress 0.1509729
## Run 17 stress 0.1023717
## Run 18 stress 0.1968676
## Run 19 stress 0.09967979
## ... Procrustes: rmse 0.000003876331 max resid 0.0000177978
## ... Similar to previous best
## Run 20 stress 0.1509889
## *** Best solution repeated 8 times
var_stress_nmds_aitch_fc <- round(nmds_source_aitch_fc$stress, 5)
var_stress_nmds_aitch_fc## [1] 0.09968
scores_source_aitch_fc= nmds_source_aitch_fc %>% vegan::scores()
meta_full_core$Depth <- factor(meta_full_core$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
nmds_plot_aitch_fc<- ggplot() +
geom_point(data=data.frame(scores_source_aitch_fc) %>%
rownames_to_column(var = "Sample")%>%
left_join(meta_full_core, by = "Sample"),
aes(x=NMDS1, y=NMDS2, color = Depth, shape = Site), size=4) +
scale_fill_lancet()+
scale_color_lancet()+
scale_shape_manual(values = c(15, 16, 17, 18, 23, 24, 25)) +
theme(axis.text = element_text(colour = "black", size = 12),
axis.title = element_text(colour = "black", size = 12),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right",
legend.box = "vertical")+ theme_bw() +
labs(title="NMDS, Aitchison distances Core") +
ylab("NMDS2")+xlab("NMDS1")
nmds_plot_aitch_fc <- nmds_plot_aitch_fc +
annotate("text", x = Inf, y = -Inf, label = paste("Stress:", var_stress_nmds_aitch_fc),
hjust = 3.1, vjust = -1.9, size = 3)
nmds_plot_aitch_fc#Combine plot
nmds_aitch_full_core <- plot_grid(nmds_plot_aitch, nmds_plot_aitch_fc, labels = c("A", "B"), ncol = 2)
nmds_aitch_full_corepseudocount <- 1
otu_table_full_core_pseudo <- otu_table_full_core + pseudocount
#Get effective reads
Effective_reads_core <- colSums(otu_table_full_core[-1, ])
meta_full_core$Effective_reads <- Effective_reads_coremeta_core_pseudo = data.frame(t(otu_table_full_core_pseudo), check.names = F) %>%
rownames_to_column(var = "Sample") %>% inner_join(meta_full_core) ## Joining with `by = join_by(Sample)`
perm_core = adonis2(t(otu_table_full_core_pseudo) ~ log(Effective_reads)+Depth_cm*Site,
data=meta_core_pseudo, method = "aitchison",
permutations = 999, strata = meta_core_pseudo$ID)
perm_coreMicrobiome
Based here
Calculate core between Depths
# count number of samples in each group
table(meta(ps)$Depth_cm, useNA = "always")##
## 0-15 16-30 31-45 50-75 <NA>
## 13 13 13 11 0
#List of Depths
Depths <- unique(as.character(meta(ps_full_rel_f)$Depth_cm))
print(Depths)## [1] "0-15" "16-30" "31-45" "50-75"
#Write a for loop to go through each of the Depths one by one and combine identified core taxa into a list.
list_core_full <- list()# an empty object to store information
for (Depth_cm in Depths){
#print(paste0("Identifying Core Taxa for ", Depth))
ps_sub_full <- subset_samples(ps_full_rel_f, Depth_cm == Depth_cm )
core_m_full <- core_members(ps_sub_full,
detection = 0, # 100 % samples
prevalence = 0.5) # 50% prevalence
print(paste0("No. of core taxa in ", Depth_cm, " : ", length(core_m_full)))
# add to a list core taxa for each group.
list_core_full[[Depth_cm]] <- core_m_full
#create core depth variable
variable_name_full <- paste0("core_full_", gsub(" ", "_", Depth_cm))
assign(variable_name_full, core_m_full, envir = .GlobalEnv)
#write all cores
file_name_full <- paste0("../results/05.diversity/",variable_name_full,".csv")
write.csv(core_m_full, file = file_name_full, row.names = FALSE)
}## [1] "No. of core taxa in 0-15 : 2048"
## [1] "No. of core taxa in 16-30 : 2048"
## [1] "No. of core taxa in 31-45 : 2048"
## [1] "No. of core taxa in 50-75 : 2048"
Plot venn diagram
library(VennDiagram)## Loading required package: grid
## Loading required package: futile.logger
##
## Attaching package: 'VennDiagram'
## The following object is masked from 'package:ggpubr':
##
## rotate
# Extract lancet Depth colors
depth_colors <- scale_fill_lancet()$palette(4)
# plot
venn_plot_full <- venn.diagram(
x = list_core_full,
filename = NULL,
#filename = "../results/05.diversity/Venn_full.png",
category.names = c("0-15", "16-30", "31-45", "50-75"),
output = TRUE,
#imagetype = "png",
height = 2000,
width = 2000,
resolution = 300,
compression = "lzw",
#units = "px",
lwd = 2,
lty = 'blank',
fontfamily ="sans",
fill = depth_colors,
cat.fontfamily = "sans",
cat.cex = 1.85,
cat.default.pos = "outer")
# Show
grid.draw(venn_plot_full)library(UpSetR)##
## Attaching package: 'UpSetR'
## The following object is masked from 'package:lattice':
##
## histogram
# Convert UpSetR fromList
full_core_data <- fromList(list_core_full)
# UpSet plot
upset_full_depth <- upset(full_core_data, nsets = 4, decreasing = T, order.by = "freq", keep.order = TRUE, main.bar.color = "black", sets.bar.color = depth_colors)
upset_full_depthGet core identities
#List of Depths
sites <- unique(as.character(meta(ps_full_rel_f)$Site))
print(sites)## [1] "Site_1" "Site_2" "Site_3" "Site_4" "Site_5" "Site_6" "Site_7"
#Write a for loop to go through each of the Depths one by one and combine identified core taxa into a list.
list_core_sites_full <- list()# an empty object to store information
for (site in sites){
#print(paste0("Identifying Core Taxa for ", Depth))
ps_sub_full <- subset_samples(ps_full_rel_f, Site == site)
core_m_full <- core_members(ps_sub_full,
detection = 0, # 100 % samples
prevalence = 0.5) # 50% prevalence
print(paste0("No. of core taxa in ", site, " : ", length(core_m_full)))
# add to a list core taxa for each group.
list_core_sites_full[[site]] <- core_m_full
#create core depth variable
variable_name_full <- paste0("core_full_", gsub(" ", "_", site))
assign(variable_name_full, core_m_full, envir = .GlobalEnv)
#write all cores
file_name_full <- paste0("../results/05.diversity/",variable_name_full,".csv")
write.csv(core_m_full, file = file_name_full, row.names = FALSE)
}## [1] "No. of core taxa in Site_1 : 2403"
## [1] "No. of core taxa in Site_2 : 3356"
## [1] "No. of core taxa in Site_3 : 2192"
## [1] "No. of core taxa in Site_4 : 3238"
## [1] "No. of core taxa in Site_5 : 2900"
## [1] "No. of core taxa in Site_6 : 2158"
## [1] "No. of core taxa in Site_7 : 2809"
Upset diagram
library(UpSetR)
sites_colors <- scale_fill_d3()$palette(7)
# Convert UpSetR fromList
input_data <- fromList(list_core_sites_full)
# UpSet plot
upset_full_sites <- upset(input_data, nsets = 7, decreasing = T, order.by = "freq", keep.order = TRUE, main.bar.color = "black", sets.bar.color = sites_colors)
upset_full_sites